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PREFACE 


This  report  presents  a user's  manual  of  operating  instructions 
for  the  "RTRAJ"  computer  program,  which  has  been  under  development  and 
in  use  at  The  Rand  Corporation  for  several  years.  (RTRAJ  takes  s 
name  from  "Rand  trajectory  program.")  The  function  of  the  program  is 
to  provide  an  estimate  of  the  accuracy  wich  which  the  motion  of  a 
vehicle  can  be  determined  on  the  basis  of  either  inertial  guidance 
measurements  or  external  tracking,  or  both. 

The  RTRAJ  program  was  developed  under  the  Future  Strategic  Aero- 
space Force  Requirements  task  of  Project  AIR  FORCE  (formerly  Project 
RAND).  It  should  be  useful  to  analysts  in  the  fields  of  ballistic 
missile  accuracy,  range  instrumentation,  and  navigational  satellites. 
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SUMMARY 


This  report  presents  the  analytical  background  and  the  operating 
instructions  ior  a computer  program,  RTRA.J,  developed  at  The  Rand  Cor- 
poration. The  purpose  of  the  program  is  to  provide  a best  estimate  of 
the  performance  of  a moving  vehicle  as  described  by  the  estimated  de- 
viation from  its  nominal  performance,  as  well  as  the  precision  with 
which  this  deviation  is  known.  The  modeled  parameters  whose  deviations 
characterize  the  vehicle  performance  include  the  vehicle  position  and 
velocity,  and  may  include  the  error  sensitivity  coefficients  associated 
with  an  inertial  measuring  unit  (1MU).  Performance  is  determined  by 
measurements  with  an  onboard  IMU  or  external  tracke.s,  or  both.  In 
the  absence  of  any  external  tracking  measurements,  the  position  and 
velocity  indicated  by  the  IMU  and  the  nominal  values  of  the  IMU  param- 
eters are  the  best  -\nd  only  estimates  of  vehicle  performance.  The  pre- 
cision of  the  position  and  velocity  estimates  is  correlated  with  the 
a priori  standard  deviations  assigned  to  the  IMU  parameters  and  varies 
as  the  vehicle  moves  along  its  path,  whereas  the  standard  deviations  of 
the  IMU  parameters  themselves  are  invariant. 

If  external  tracking  measurements  of  known  precision  are  made  in 
conjunction  witn  the  IMU,  this  additional  information  can  be  used  to 
modify  the  performance  estimate  provided  by  the  IMU  by  means  of  a 
sequential  Kalman  filter,  which  produces  an  improved  estimate  of  the 
vehicle  position  and  velocity  as  well  as  of  the  IMU  parameters.  At  the 
same  time,  the  precision  of  these  estimates  is  also  redetermined.  In 
the  event  that  no  IKU  is  used,  the  external  measurements  are  used  to 
update  the  a priori  estimates  of  performance  and  precision.  In  the 
absence  of  any  unmodeled  errors,  the  best  estimate  of  the  vehicle  per- 
formance should  approach  its  actual  performance  with  increasing  preci- 
sion as  the  number  of  external  measurements  increases.  It  is  possible, 
however,  to  introduce  unmodeled  errors  in  the  form  of  unknown  tracker 
position  and  velocity  biases  as  well  as  biases  in  the  measured  quanti- 
ties. With  such  unmodeled  errors,  the  estimated  vehicle  performance 
deviates  from  the  actual  performance. 
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The  trackers  used  to  make  the  external  measurements  can  be  lo- 
cated on  the  ground,  on  aircraft,  on  satellites,  or  on  the  vehicle  (as 
an  altimeter),  and  can  measure  any  or  all  of  the  following  quantities: 

! 

range,  range  rate,  azimuth  angle,  elevation  angle,  hyperbolic  range,  i 

and  hyperbolic  range  rate. 

H 

1 

1 

3 


i 

\ 

\ 


i 

3 

1 


-vii- 


ACKNOWLEDGMENTS 


Several  Rand  staff  members  have  been  instrumental  in  the  develop- 
ment of  the  RTRAJ  program.  They  include  R.  L.  Mobley,  who  wrote  the 
original  version;  R.  P.  Castro,  who  incorporated  the  inertial  guidance 
computations;  and  Eric  Olson,  who  added  the  hyperbolic  range  and  range 
rate  measurements.  The  author  would  like  to  acknowledge  their  assis- 
tance in  translating  the  program  from  the  original  FORTRAN  into  the 
language  of  the  present  report. 


CONTENTS 


PREFACE  iii 

SUMMARY  v 

ACKNOWLEDGMENTS  vii 

FIGURES  xi 

SYMBOLS  xiii 


Section 

I.  INTRODUCTION  1 

II.  ANALYSIS  5 

Vehicle  Motion  5 

Inertial  Measuring  Unit  Error  Analysis  5 

The  State  Vector  - 9 

The  State  Transits  i Matrix  (Vehicle  Model)  10 

Covariance  Matrix  (Propagation  Through  Time)  12 

Modeling  the  External  Measurement  Variables  13 

Estimation  of  State  Vector  and  Covariance  Matrix  ..  14 

III.  PROGRAM  OPTIONS  16 

Mode  of  Operation 16 

Flight  Path  Characteristics  19 

Tracker  Characteristics  20 

Measureme it  Characteristics  21 

IV.  MEASUREMENT  PROCEDURE  23 

Tracker  Selection  23 

Measurement  Types  25 

V.  RTRAJ  PROGRAM  OPERATION 36 

Input  Deck  „ 36 

Output  Format  53 

Appendix 

COORDINATE  SYSTEMS  61 


#PRSC«Dm  -AfiS.'HJJCC.NOT  flV&D 


FIGURES 


Earth-centered  rotating  coordinate  system  62 

Launch-centered  rotating  coordinate  system  63 

In-plane  rotating  coordinate  system  64 

Accelerometer  coordinate  system  67 

Tracker  coordinate  system  68 


\ 5***52®  *■«* . SUMUM. 


-xii  1- 


SYMBOLS 
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[Ap] 
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Matrix  transformation  between  [£]  and  [H]. 

Matrix  transformation  from  earth-centered  inertial 
(ECI)  to  in-plane  inertial  (IPX)  coordinates. 

Matrix  transformation  from  ECI  to  in-plane  rotating 
(IPR)  coordinates. 

Components  of  Aa  along  the  input,  cross,  and  normal 
axes  of  the  ith  accelerometer. 

Components  of  Aa  along  the  input,  output,  and  spin 
axes  of  the  ith  gyroscope. 

Matrix  transformation  from  ECI  to  launch-centered 
inertial  (LCI)  coordinates. 

Matrix  transformation  from  ECI  to  launch-centered 
rotating  (LCR)  coordinates. 

Matrix  transformation  from  ECI  to  platform  coordinates. 

Matrix  transformation  from  platform  coordinates  to 
misaligned  platform  coordinates. 

Matrix  transformation  from  earth-centered  rotating 
(ECR)  to  the  ith  tracker  coordinates. 

Nongravitational  acceleration  of  the  missile. 

Matrix  transformation  from  platform  coordinates  to 
the  ith  accelerometer  coordinates. 

Matrix  transformation  from  platform  coordinates  to 
the  ith  gyroscope  coordinates. 

Vector  of  the  errors  in  the  sensed  accelerations  of 
the  three  accelerometers. 

Vector  of  the  acceleration  error  due  to  the  k th  drift 
rate  component  of  the  ith  gyroscope. 

Vector  of  the  acceleration  errors  due  to  platform 
misalignment . 

Total  vector  error  in  the  sensed  acceleration. 

Component  of  AA^  for  the  ith  accelerometer. 

Matrix  of  the  partial  derivatives  or  the  measured 
quantity  with  respect  to  the  missile's  coordinates 
relative  to  the  tracker. 

Coefficients  In  the  expressions  for  the  standard  devi- 
ations of  the  external  measurements,  where  i runs  from 
1 to  24. 
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Components  ot  gravity  along  the  input,  cross,  and 
normal  axes  of  the  ith  accelerometer. 


Si’Si’Si 


Components  of  gravity  along  the  input,  output,  and 
spin  axes  of  the  ith  gyroscope. 


(H)  Parameter  transformation  matrix. 
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State  covariance  matrix. 
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determined  by  the  IMU. 
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Random  variation  in  the  range  measurement. 

Random  variation  in  the  range  rate  measurement. 

Bias  in  the  range  measurement . 

Bias  in  the  range  rate  measurement . 
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Best  estimate  of  the  sta^e  vector. 
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Standard  deviation  of  the  azimuth  angle  measurement. 
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I.  INTRODUCTION 


The  RTRAJ  program  described  in  this  report  is  designed  to  produce 
a best  estimate  of  the  performance  of  a moving  vehicle,  based  on  measure- 
ments made  by  an  onboard  inertial  measuring  unit  (IMU)  or  external 
trackers,  or  both.  The  performance  of  a vehicle  can  be  specified  by  a 
set  of  state  variables  which  includes  the  vehicle  position  and  velocity, 
and  may  include  the  error  sensitivity  coefficients  of  an  onboard  IMU. 

For  the  purposes  of  the  analysis,  it  is  convenient  to  define  these  state 
variables  by  their  differential  deviations  from  a nominal  flight  path 
and  nominal  IMU  characteristics.  It  is  the  purpose  of  the  program  to 
make  a best  estimate  of  these  differential  variations,  which  constitute 
the  state  vector  of  the  system.  In  addition,  the  precision  of  the  esti- 
mate of  the  elements  of  the  state  vector  is  also  determined,  in  the  form 
of  a covariance  matrix  whose  diagonal  elements  are  the  variances  of  the 
state  variables  and  whose  off-diagonal  terms  are  a measure  of  the  co- 
variance  between  state  variables.  In  such  a system,  the  position  and 
velocity  indicated  by  the  IMU  is,  in  the  absence  of  any  external  mea- 
surements, the  best  and  only  estimate  of  the  actual  position  and  velocity 
of  the  vehicle.  Under  these  conditions,  the  estimates  of  the  deviations 
of  the  state  variables  from  their  nominal  values  f*re  identically  zero. 
However,  the  accelerations  sensed  by  the  IMU  are  functions  of  the  state 
variables  associated  with  the  IMU,  with  the  result  that  a priori  stand- 
ard deviations  in  these  IMU  parameters  are  propagated  into  standard  de- 
viations in  position  and  velocity  as  time  advances,  with  a corresponding 
modification  of  the  position  and  velocity  elements  of  the  covariance 
matrix.  At  the  same  time,  the  standard  deviations  of  the  IMU  parameters 
remain  constant. 

If,  in  addition  to  the  IMU  output,  measurements  are  available  from 
external  trackers  in  the  form  of  range,  range  rate,  azimuth,  or  eleva- 
tion angle  determinations,  these  data  can  be  used  to  improve  the  estimates 
of  the  state  variables  and  their  standard  deviations.  This  updating  of 
the  estimates  is  accomplished  by  means  of  a sequential  Kalman  filter, 
which  combines  the  current  best  estimate  of  the  state  variables  with  a 
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given  measured  value,  the  standard  deviation  of  the  measurement,  and 
the  sensitivity  of  the  state  variables  to  this  particular  measurement 
to  provide  a new  best  estimate  of  the  state  variables.  At  the  same 
time,  these  measurement  data  are  also  used  to  update  the  standard 
deviations  of  the  estimates  of  the  state  variables.  This  procedure 
enables  the  user  to  estimate  the  deviation  of  the  system  from  its  nomi- 
nal performance  and  to  improve  his  confidence  in  this  estimate.  It 
should  be  noted  that  use  of  this  Kalman  filter  assumes  that  the  indi- 
vidual measurements  are  processed  sequentially  and  that  they  are  un- 
correlated. In  addition,  if  unraodeled  errors  are  present,  such  as 
unknown  biases  in  tracker  positions  or  in  the  measurements,  the  result- 
ing estimate  of  the  vehicle  performance  will  deviate  from  the  actual 
performance.  Although  such  biases  could  be  modeled  as  additional  state 
variables  to  be  estimated,  this  option  is  not  available  in  the  present 
version  of  the  program. 

As  an  example  of  the  type  of  problem  to  which  the  program  is  ap- 
plicable, consider  the  flight  of  a ballistic  missile  with  a period  of 
powered  flight  followed  by  an  interval  of  free  fall  to  impact.  It  is 
assumed  that  vehicle  position  and  velocity  are  computed  by  an  onboard 
IMU  and  that  this  information  is  telemetered  to  ground  observers.  Dur- 
ing powered  flight,  the  motion  of  the  missile  is  also  observed  by  ground 
trackers  which  each  measure  the  instantaneous  range  to  the  missile  at 
specified  measurement  intervals.  Initially,  with  the  missile  on  the 
launch  pad,  there  is  no  deviation  in  either  position  or  velocity  from 
that  of  the  nominal  flight  path.  However,  the  error  sensitivity  param- 
eters of  the  particular  IMU  may  deviate  from  their  nominal  values. 
Similarly,  the  a priori  standard  deviations  of  the  position  and  velocity 
are  initially  zero,  whereas  those  of  the  IMU  parameters  are  not.  When 
the  missile  is  launched,  it  flies  a path  such  that  the  accelerations 
sensed  by  the  IMU  are  those  specified  for  the  nominal  trajectory,  but, 
since  the  IMU  is  in  error,  the  actual  missile  deviates  from  the  nominal 
trajectory.  As  indicated  above,  if  the  IMU  were  the  only  source  of 
information,  this  deviation  from  the  nominal  flight  path  could  not  be 
detected.  However,,  when  the  tracking  information  is  incorporated  into 
the  estimation  procedure,  it  is  possible  not  only  to  obtain  an  estimate 


of  the  deviations  in  position  and  velocity  from  the  nominal  trajectory 
but  also  of  the  nonstandard  values  of  the  IMU  parameters.  In  addition, 
an  improved  estimate  of  the  standard  deviations  of  these  estimates  of 
the  state  variables  is  achieved.  If  no  further  measurements  are  made 
with  either  the  IMU  or  the  external  trackers,  the  state  variables  and 
the  associated  covariance  matrix  can  be  propagated  from  their  values 
at  cutoff  along  the  trajectory  to  impact.  The  estimated  and  actual  dis- 
placement of  the  impact  point  from  that  of  the  nominal  trajectory  can 
then  be  determined,  as  well  as  the  dimensions  and  orientation  of  the 
error  ellipse  about  that  point. 

Another  example  is  tracking  a ballistic  missile  during  that  por- 
tion of  its  flight  from  cutoff  to  impact.  In  this  case,  there  is  no 
IMU  operating,  but  a set  of  ground  trackers  is  measuring  the  range  rate 
of  the  missile  relative  to  each  trac’  er.  It  is  now  necessary  to  specify 
a nominal  initial  position  and  velocity  of  the  missile  at  cutoff,  which 
determines  the  nominal  flight  path.  In  addition,  the  initial  deviations 
of  the  actual  missile  position  and  velocity  from  these  nominal  condi- 
tions, as  well  as  their  a priori  standard  deviations,  can  be  introduced. 
As  in  the  previous  case,  the  actual  deviations  from  the  nominal  flight 
path  are  determined,  as  well  as  the  best  estimate  of  these  deviations, 
based  on  combining  the  a priori  estimates  with  the  tracker  data  by  means 
of  the  Kalman  filter.  In  this  «ay,  the  user  obtains  an  estimate  of  the 
location  of  the  actual  impact  as  well  as  the  confidence  in  this  estimate. 

A final  example  is  the  application  to  a navigational  satellite 
system  in  which  trackers  are  mounted  on  satellites  whose  ephemerides 
are  known.  As  in  the  previous  example,  the  user  vehicle,  which  could 
be  a ballistic  missile,  another  satellite,  an  aircraft,  or  a terres- 
trial vehicle,  is  observed  by  the  trackers  and  its  estimated  and  actual 
deviation  from  its  specified  nominal  path  are  determined,  together  with 
the  precision  of  this  estimate. 

This  report  is  intended  as  a user's  manual  for  the  operation  of 
the  RTRAJ  program  and  not  as  a programmer's  manual,  and  thus  it  does 
not  contain  the  program  listing  or  flow  charts.  The  organization  of 
the  report  is  as  follows:  Section  II  presents  the  analysis  necessary 

'.o  describe  the  vehicle  motion,  the  inertial  measuring  unit,  the  propa- 
gation of  the  state  vector  and  the  state  covariance  matrix,  and,  finally, 


the  operation  of  the  Kalman  filter  to  update  the  performance  estimate 
on  the  basis  of  the  external  measurements.  Section  III  describes  the 
various  options  available  in  the  program.  Section  IV  delineates  the 
various  types  of  tracking  measurements  available  and  how  they  are  pro- 
cessed to  provide  input  for  the  Kalman  filter.  Section  V discusses  the 
operation  of  the  program,  including  both  the  configuration  of  the  input 
deck  and  the  various  output  formats  which  can  be  invoked.  Finally,  the 
appendix  outlines  the  various  coordinate  systems  used  in  the  program 
and  the  transformation  matrices  relating  them. 

The  program  is  written  in  FORTRAN  IV  and  has  been  implemented  on 
an  IBM  370/158  computer. 
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II.  ANALYSIS 


VEHICLE  MOTION 

The  nominal  motion  of  a vehicle  in  the  earth's  gravitational  field, 
is  defined  by  the  equation 


R 

m 


m 


(1) 


where  R is  the  vehicle  position  vector  relative  to  the  earth's  center, 

U is  the  earth's  gravitational  constant,  and  A is  the  nongravitational 

& 

acceleration  including  thrust  and  aerodynamic  forces. 

The  motion  of  the  vehicle,  as  indicated  by  the  IMU,  is  defined  by 
the  equation 


JL*  MRr,  _ _ 

^ . + * + A* 

Rp 


(2) 


where  Rp  is  the  vehicle  position  vector  as  indicated  by  the  IMU,  and 

AA  is  the  acceleration  error  resulting  from  uncertainties  in  the  ac- 
0 

celerometer  parameters  as  well  as  platform  misalignment  and  drift. 

INERTIAL  MEASURING  UNIT  ERROR  ANALYSIS 

The  motion  of  the  vehicle  is  indicated  by  an  onboard  inertial 
measuring  unit,  which  consists  of  three  accelerometers  fixed  to  a 
stable  platform  with  a 56-deg  angle  between  the  sensitive  axes  of  each 
pair  of  accelerometers.  The  stabilization  of  the  platform  is  accom- 
plished by  three  single-degree-of-f reedom  gyroscopes  mounted  with 
mutually  perpendicular  input  axes.  This  configuration  is  used  unless 
other  instrument  orientation  angles  are  specified  in  the  input  deck. 

It  is  assumed  that  the  input  acceleration  is  expressed  in  platform 
coordinates  even  though  these  coordinates  may  be  misaligned  relative 

to  their  idealized  position.  This  input  acceleration  is  represented 

-(P) 

as  a column  vector,  A 

a 


Accelerometer  Errors 


It  is  convenient  for  discussion  to  project  the  applied  accelera- 
tion into  accelerometer  axes  by  the  following  transformation: 
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where  A is  the  acceleration  in  the  instrument  coordinates  of  the 
a 

1th  accelerometer  and  A^,  A^,  and  A^  are  the  components  along  the 
input,  cross,  and  normal  axes,  respectively,  of  the  1th  accelerometer. 

There  are  six  sources  of  error  associated  with  each  accelerometer 
the  error  in  the  sensed  acceleration  for  the  1th  accelerometer  is 
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where  the  K's  are  the  error  sensitivity  coefficients,  and  the  y's  are 
given  as 
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where  GT,,  G_ . , and  G„.  are  the  components  of  gravity  along  the  I,  C, 
It  Cl  Nl 

and  N axes,  respectively,  just  prior  to  launch.  The  G's  are  only  used 
when  gravity  compensation  is  used;  otherwise,  they  are  set  to  zero. 


The  errors  in  the  sensed  accelerations  for  the  three  accelerometers 
can  be  represented  as  a column  vector  of  the  form 


(5) 


which  can  be  transformed  to  platform  coordinates  by  the  relation 


[R]“1A^1) 


(6) 


where  R is  the  mrtrix  transformation  from  platform  coordinates  to  the 
three  accelerometer  input  axes.  Since  these  axes  are  not  orthogonal, 
it  is  necessary  to  use  the  inverse  of  R to  accomplish  the  transforma- 
tion in  Eq.  (6). 


Platform  Misalignment  Errors 

If  the  inertial  platform  is  initially  misaligned  by  small  angular 
rotations  of  K^,  1^,  and  about  the  y,  z,  and  x platform  axes,  respec- 
tively, the  resulting  acceleration  errors  introduced  along  the  unper- 
turbed platform  axes  are  given  by  the  relation 


(7) 


Gyroscope  Drift  Errors 

In  addition  to  the  initial  platform  misalignment,  the  drift  of 
the  three  stabilizing  gyroscopes  may  result  in  further  misalignment  as 
time  goes  on.  To  determine  the  drift  rates  of  these  gyroscopes,  it  is 
again  convenient  to  project  the  applied  acceleration  along  the  instru- 
ment axes  by  the  relation 
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where  \J^,  A^,  and  are  the  acceleration  components  along  the  in- 
put, output,  and  spin  axes,  respectively,  ot  the  i t.h  gyroscope. 

There  are  eight  sources  of  gyroscope  drift;  the  drift  rate  of  the 
1th  gyroscope  due  to  the  k t.h  source  of  drift  is  given  by 

4>ki  = K3(k-2)+iyki  Ot  = 1 through  16)  , (9) 


where,  as  in  the  case  of  the  accelerometers,  the  K's  are  the  error 
sensitivity  coefficients,  and  the  p's  are  given  by  the  following 
relations: 
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where  G^,  G^,  and  G are  the  components  of  gravity  along  the  input, 
output,  and  spin  axes,  respectively,  immediately  prior  to  launch.  As 
in  the  case  of  the  accelerometers,  these  G terms  are  omitted  if  gravity 
compensation  is  not  used.  The  last  term  in  Eq.  (9)  is  a time-dependent 
drift  rate  and  is  proportional  to  the  elapsed  time  since  launch,  T - T . 
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As  a result  of  these  drift  rates,  small  angular  rotations  of  the 
platform  develop  about  the  three  mutually  perpendicular  gyroscope 
input  axes  and  are  given  by 


(10) 


The  resulting  acceleration  error ( in  unperturbed  platform  coordinates 
due  to  the  kth  drift  rate  of  the  ith  gyroscope  is  given  by 
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Combined  Acceleration  Errors 

All  of  the  acceleration  errors  discussed  above  make  up  the  quan- 
tity Aa  1..  Eq.  (2).  Since  this  equation  is  expressed  in  ECI  coordi- 
nates,  the  value  of  the  acceleration  error  term  is  given  by 
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This  expression  involves  a total  of  45  error  sensitivity  parameters, 
represented  by  the  K's  which,  for  an  ideal  system,  would  all  be  zero 
with  the  exception  of  the  three  scale  factor  parameters,  which 
should  be  unity. 


THE  STATE  VECTOR 

The  state  vector  dX  that  describes  this  system  has  51  components, 
including  the  small  variations  from  nominal  of  the  three  missile  posi- 
tion components,  the  three  missile  velocity  components,  and  the  45 
1MU  parameters  described  above.  Thus,  given  the  nominal  flight  path 
and  nominal  IMU  parameters,  the  system  performance  can  be  determined 
by  solving  for  the  state  vector  as  a function  of  time,  as  dX(t). 


The  state  vector  at  time  t is  related  to  its  value  at  time 
through  the  state  transition  matrix  ( t , t q)  by  the  relation 

dX(t)  - [i(t,tQ) ] dX(t0)  . (13 

The  differential  equation  for  i(t,t.)  is  obtained  by  first  different! 
ating  Eq.  (13)  with  respect  to  time  to  give 

dX(t)  - [*(t,t0)]  dX(tQ)  . (14 

Elimination  of  dX(t^)  between  Er  (13)  and  (14)  results  in 

dX(t)  - [*(t,t0)]  [*(t,t0)]"1  dX(t)  . (15 


which  is  the  desired  differential  equation  for  ♦(t.t^),  and  where 


The  matrix  [A]  is  a square  matrix  which  can  be  partitioned  as  follows 
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The  partitioned  elements  can  be  evaluated  as 
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“3K  1 

— - - [0]  (45  x 3) 
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= [0]  (45  x 3) 

_3Rp_ 
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The  final  (3  x 45)  element,  vR^/SK^,  is  evaluated  by  differenti- 
ating Eq.  (2)  with  respect  to  che  45  IMU  error  sensitivity  parameters 

in  the  expression  defining  AA  . 

3 

Since  the  lower  45  rows  of  the  A matrix  are  identically  zero,  it 
follows  from  Eq.  (17)  that  the  corresponding  45  rows  in  the  $ matrix 
are  also  zero.  Since  the  $ matrix  is  initially  an  identity  matrix 
when  t equals  tQ,  its  lower  45  rows  remain  unchanged  with  time  because 
the  differential  values  of  the  IMU  parameters  remain  fixed.  Thus,  to 
determine  the  new  value  of  the  state  transition  matrix,  it  is  only 
necessary  to  determine  the  elements  in  the  upper  six  rows  of  4>,  which 
involves  the  integration  of  306  first-order  differential  equations  in 
the  form  of  Eq.  (17).  In  addition,  there  are  24  equations  of  the  form 
of  Eq.  (9)  which  must  be  integrated  to  determine  the  drift  angles  asso- 
ciated with  the  platform,  and,  finally,  six  first-order  equations 
resulting  from  q.  (1)  to  give  the  three  components  of  missile  position 
and  velocity.  Thus,  a total  of  336  first-order  differential  equations 
must  be  integrated  to  determine  the  nominal  missile  trajectory  and  the 
state  transition  matr*:  . The  state  transition  matrix  can  then  be  used 
in  Eq.  (13)  to  propagate  the  state  vector. 


COVARIANCE  MATRIX  (PROPAGATION  THROUGH  TIME) 

The  covariance  matrix  associated  with  the  state  vector  is  a sym- 
metrical 51  x 51  matrix  where  the  diagonal  elements  equal  the  squares 
of  the  standard  deviations  or  variances  of  the  state  variables,  and 
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where  the  off-diagonal  elements  are  the  covariances  indicating  the 
correlation  between  variables.  Initially,  the  first  six  diagonal  ele- 
ments representing  the  variances  of.  the  three  position  and  the  three 
velocity  coordinates  of  the  vehicle  are  set  equal  to  zero,  since  then, 
is  no  uncertainty  in  position  or  velocity  prior  to  launch.  The  other 
45  diagonal  elements  are  set  equal  to  the  variances  of  the  IMU  param- 
eters, and  the  off-diagonal  covariance,  terms  are  all  zero.  As  the 
flight  progresses,  the  covariance  matrix  can  be  propagated  by  means  of 
the  state  transition  matrix  by  the  relation 

[Q(tk)l  - [*(tk,tJ)]  [QO^)]  [*(tk,tj)lT  , (20) 

which  computes  the  covariance  matrix  Q(tk>  at  time  tk  from  its  previous 
value  Q(tj).  Because  of  the  form  of  the  state  transition  matrix  <J>, 
the  propagated  form  of  the  covariance  matrix  Q(tk)  can  develop  nonzero 
values  in  the  elements  in  the  first  six  rows  or  first  si.  columns,  or 
both.  Thus,  the  variances  of  the  IMU  parameters  result  in  uncertain- 
ties in  the  measured  values  of  the  missile's  position  and  velocity. 
Although  this  process  results  in  correlation  between  the  IMU  parameters 
and  both  position  and  velocity,  there  is  no  correlation  between  the 
various  IMU  parameters,  as  evidenced  by  the  fact  that  the  lower  right 
hand  45  * 45  element  of  the  covariance  matrix  remains  unchanged — having 
all  elements  zero  except  those  on  the  diagonal,  which  remain  equal  to 
the  fixed  variances  of  the  IMU  parameters. 

MODELING  THE  EXTERNAL  MEASUREMENT  VARIABLES 

The  foregoing  analysis  has  considered  only  the  effect  of  the  vari- 
ations of  the  IMU  parameters  on  the  accuracy  of  the  determination  of 
the  missile's  position  and  velocity.  If,  however,  external  measure- 
ments of  the  missile's  motion  are  made  with  known  precision,  during 
its  flight,  these  measurements  caw  be  combined  with  those  of  the  guid- 
ance system  to  provide  an  improved  estimate,  not  only  of  the  errors  in 
missile  position  and  velocity  but  also  of  the  errors  in  the  IMU  param- 
eters. It  is  first  necessary,  however,  to  determine  the  most  likely 


value  of  the  measurement  errors  based  on  the  current  best  estimate  of 
the  state  vector.  This  is  given  by  the  relation 


dZg  = [H]  dX  , (21) 

where  dZg  is  the  most  likely  measurement  error,  and  [H]  is  the  param- 
eter transformation  matrix  with  a dimension  of  1 x 51  and  is  repre- 
sented as 
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Since  the  external  measurement  is  independent  of  the  IMU  parameters, 
the  last  45  elements  of  L H ] are  zero.  The  expressions  for  the  remain- 
ing six  elements  of  the  H matrix  for  various  types  of  measurements  are 
derived  in  Sec.  IV  of  this  report. 

ESTIMATION  OF  STATE  VECTOR  AND  COVARIANCE  MATRIX 

By  means  of  a Kalman  filter,  the  external  measurements  can  be  com- 
bined with  observations  of  the  inertial  guidance  system  to  provide  an 
improved  estimate  of  the  elements  of  the  state  vector  as  well  as 
of  the  state  covariance  matrix.  Initially,  at  time  t^,  the  state 
transition  matrix  is  set  equal  to  an  identity  matrix  and  the  integra- 
tion of  the  equations  of  motion  proceeds  to  time  t^,  when  the  first 
measurement  is  me  e.  At  this  point,  the  integration  has  generated  the 
state  transition  matrix  $(t^,tQ),  which  is  now  used  to  propagate  the 
state  vector  and  the  state  covariance  to  time  t^  by  means  of  Eqs.  (13) 
and  (20).  These  propagated  matrices  are  represented  as  dX(t^/0)  and 
Q (t^/0) , where  the  zero  in  the  argument  means  that  no  measurement  data 
have  yet  been  incorporated. 

These  matrices  can  now  be  updated  by  means  of  the  Kalman  filter 
relations  as  follows:* 

*R.  H.  Battih,  Astronautiaal  Guidance , McGraw-Hill  Book  Company, 
Inc. , New  York,  1964. 
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(Q(t1/0)]  [Hj1  TdZ.  - (H  ] [dX(t,/0)]l 

dX(t./l)  = dX(t  /O)  + - - k— ^ 1 (23) 

IHj]  [Q(t1/0)1  [H^1  + 01 

lQ(t./0)]  [HJT  [H.  1 tQ(t , /O)  ] 

[Q(t./l)l  = [Q(t. /O) ] ± i — — j (24) 

1 IH^  tQ(t1/0)  ] + oj 

where  dZ^  is  the  difference  between  the  actual  measurement  and  that 
of  the  nominal  reference  trajectory,  [H^]  is  the  H matrix  row  corre- 
sponding to  measurement  number  one,  and  is  the  standard  deviation 
of  measurement  number  one. 

If  a total  of  n measurements  are  made  at  time  t^,  Eqs.  (23)  and 
(24)  are  applied  n times,  after  which  the  best  estimates  of  the  state 
vector  and  the  state  covariance  are  dX(t^/n)  and  Q(t^/n).  After  all 
of  the  measurements  have  been  processed,  the  state  transition  matrix 
^>^tl’tl^  *S  Set  e<*ua*  to  an  ident*ty  matrix  and  the  above  procedure 
is  repeated  until  missile  impact  occurs.  At  this  time,  the  first  three 
elements  of  the  state  vector,  representing  the  deviations  in  the  esti- 
mated missile  position  from  the  nominal  trajectory  impact  point,  are 
used  to  determine  the  displacement  of  the  estimated  impact  point.  In 
addition,  the  state  covariance  matrix  can  be  used  to  determine  the 
major  and  minor  axes,  as  well  as  the  orientation,  of  the  error  ellipse 
on  the  ground  about  the  estimated  Impact  point. 
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III.  PROGRAM  OPTIONS 


Since  the  computer  program  described  in  this  report  can  be  applied 
to  a wide  variety  of  problems,  there  are  a number  of  options  which  can 
and  must  be  specified  in  any  given  run.  These  options  are  described 
below;  the  mechanics  of  invoking  them  are  shown  in  Sec.  V. 

MODE  OF  OPERATION 

In  Sec.  II,  it  is  shown  that  the  nominal  performance  of  a missile 
is  specified  by  the  solution  of  Eq.  (1),  which  assumes  no  errors  in 
the  IMU  parameters.  The  nonnominal  performance  of  this  same  missile 
can  be  characterized  by  a state  vector  consisting  of  the  position  and 
velocity  errors  relative  to  the  nominal  trajectory,  as  well  as  the  45 
nonzero  values  of  the  £MU  parameters.  Before  discussing  the  modes  of 
operation,  it  is  convenient  to  consider  two  forms  of  this  state  vector. 
The  first,  dXc>  describes  the  actual  performance  of  the  missile,  since 
its  elements  are  the  actual  deviations  in  position  and  velocity  from 
the  nominal  trajectory  and  the  actual  values  of  the  errors  in  the  IMU 
parameters.  These  elements  of  dX£  can  be  selected  from  gaussian  dis- 
tributions with  appropriately  specified  standard  deviations.  In  any 
event,  these  values  are  not  available  to  the  onboard  guidance  system. 

The  second  state  vector,  dX,  represents  the  performance  of  the  missile 
as  estimated  by  the  onboard  guidance,  with  or  without  external 
measurements. 

Mode  0 (Propagation  of  Covariance) 

In  this  case,  it  is  assumed  that  the  missile  flies  the  nominal 
trajectory;  the  elements  of  dX£  are  thus  identically  zero.  Initially, 
all  of  the  elements  of  the  state  vector  dX  are  also  zero,  since  at  this 
time  there  is  no  reason  to  doubt  the  IMU  output.  In  the  covariance 
matrix,  the  only  nonzero  values  are  the  diagonal  elements  representing 
the  variances  of  the  IMU  parameters.  As  the  missile  progresses  along 
its  trajectory,  in  the  absence  of  any  external  measurement,  the  values 
of  dX  and  dXc  remain  zero,  and  the  state  covariance  matrix  is  propagated 
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in  accordance  with  Kq.  (20).  If  external  measurements  are  made,  the 
values  of  dX  .nay  still  remain  zero  unless  there  are  biases  in  the  mea- 
surements. In  the  presence  of  a measurement  bias,  the  term  dZ^  in 
Eq.  (23)  is  no  longer  zero  and  the  resulting  values  of  dX  are  likewise 
nonzero.  At  the  same  time,  the  covariance  matrix  is  modified  by  the 
measurement  process  of  Eq.  (24)  to  give  an  improved  estimate  of  the 
variances  of  the  state  variables  as  well  as  the  off-diagonal  covar- 
iances. It  should  be  emphasized  that  the  resulting  covariance  matrix 
is  independent  of  the  biases  in  the  measurement  and  that  the  nonzero 
values  of  dX  represent  an  erroneous  estimate  of  the  missile  performance, 
since  the  missile  actually  is  on  the  nominal  trajectory  (dXc  is  zero). 

This  mode  of  operation  can  also  be  used  to  estimate  the  nonnominal 
performance  of  the  missile  as  characterized  by  the  state  vector  dXc . 

To  do  this,  the  initial  position  and  velocity  errors  relative  to  the 
nominal  trajectory,  as  well  as  the  1MU  parameters,  are  selected  from 
gaussian  distributions  to  form  the  state  vector  dX( , where  the  state 
vector  dX  is  initially  zero.  Both  dX£  and  dX  are  propagated  by  means 
of  Eq.  (13)  up  to  the  time  of  the  first  external  measurements,  when 
dX  is  still  zero.  The  value  of  dX  is  then  updated  by  incorporating 
the  external  measurements  by  means  of  Eq.  (23),  after  which  dX  is  again 
propagated  to  the  time  of  the  next  measurements.  As  this  procedure  of 
propagating  and  updating  proceeds,  the  estimate  of  the  state  vector 
dX  should  approach  the  actual  state  vector  dXc , thus  giving  an  improved 
estimate  of  the  missile's  performance.  The  alternate  propagation  and 
updating  of  the  covariance  matrix  by  means  of  Eqs.  (20)  and  (24)  gives 
the  same  result  as  that  obtained  abc  e,  since  this  propagation  and 
updating  is  independei.  the  actual  state  vector  dXc> 

Mode  1 (Generation  of  Data) 

This  mode  of  operation  is  designed  to  produce  a set  of  tracking 
data  that  can  be  used  later  as  input  to  the  Mode  2 operation  of  the 
program,  in  which  an  estimate  of  the  missile  performance  is  determined. 
The  data  so  generated  are  of  the  same  sort  that  might  be  provided  by 
tracking  stations  during  an  actual  missile  firing. 
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To  generate  the  appropriate  tracking  data,  the  missile  flies  along 
the  nominal  trajectory  specified  by  Eq.  (1).  Each  of  the  external  track- 
ing stations  is  assumed  to  be  capable  of  measuring  one  or  more  of  the 
following  quantities;  the  range  to  the  missile,  its  range  rate,  its 
azimuth  angle,  or  its  elevation  angle.  At  specified  time  , for  each  of 
the  tracking  stations,  the  true  values  of  these  observables  art  deter- 
mined. These  values  are  then  modified  by  introducing  biases  in  the 
tracker  position  and  in  the  measurement,  and,  finally,  a andom  measure- 
ment error  selected  from  a gaussian  distribution  whose  standard  devia- 
tion is  that  specified  for  the  particular  measurement  process.  These 
modified  values  are  the  desired  observational  data  and  are  recorded  to- 
tether  with  the  tracker  identification  and  the  time  of  the  measurement. 

Mode  2 (Simulation) 

In  this  mode  of  operation,  the  objective  is  to  estimate  the  mis- 
sile state  vector  dX  by  combining  the  IMU  measurements  with  the  external 
tracking  data.  As  indicated  above,  these  external  tracking  data  can  be 
provided  either  by  the  Mode  1 operation  or  by  recording  the  tracking 
data  from  an  actual  missile  flight.  Initially,  the  elements  of  dX,  the 
estimated  state  vector,  are  set  equal  to  zero  since,  in  the  absence  of 
any  external  measurement,  there  is  no  reason  to  assume  that  the  vehicle 
has  deviated  from  the  nominal  trajectory.  However,  at  each  measurement 
time,  the  external  measurements  are  used  to  update  the  estimated  state 
vector  dX.  These  values  are  then  used  to  determine  the  new  values  of 
the  state  variables  by  the  relation 

X(t/n)  = X(t/m)  + dX(t/n)  , (25) 

where  the  arguments  of  X have  the  same  meaning  as  those  of  dX  in  Eq. 

(23),  and  the  m and  n are  the  number  of  measurements  before  and  after 
updating,  respectively.  Equation  (25)  results  in  a specification  of 
a new  nominal  trajectory  which  should  approach  the  true  trajectory  es 
the  number  of  external  measurements  increases.  It  should  be  noted  t.ii._ 
this  type  of  simulation,  using  Mode  1 to  generate  the  data,  is  essen- 
tially equivalent  to  the  Mode  0 operation,  with  a nonzero  value  of  dXc. 
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Thus,  the  principal  advantage  of  Mode  2 is  its  ability  to  utilize 
externally  generated  real  tracking  data.  Otherwise,  Mode  0 is  prefer- 
able and  uses  less  computer  time  per  run. 

FLIGHT  PATH  CHARACTERISTICS 


Staging 

For  computational  convenience,  the  flight  path  is  divided  into 
several  time;  intervals  or  stages.  At  the  beginning  of  each  stage,  data 
are  input  to  describe  the  vehicle  acceleration,  the  characteristics  of 
the  external  tracking  system,  and  the  format  and  frequency  of  the  print- 
out. Thus,  a given  run  might  have  a stage  of  powered  flight  in  which 
the  vehicle  acceleration  is  specified  and  during  which  no  tracking  mea- 
surements are  made.  This  could  be  followed  by  a stage  of  unpowered 
flight  in  which  the  missile  is  in  free  fall  and  during  which  external 
measurements  are  made  by  the  ground  trackers.  Finally,  this  would  be 
followed  by  another  stage  of  free  fall,  but  with  no  external  measure- 
ments, continuing  on  to  impact. 

Acceleration  Profile 

The  acceleration  of  the  missile  as  a function  of  time  can  be  input 
in  tabular  form  by  expressing  the  three  components  of  acceleration  at 
specified  times.  The  acceleration  at  any  other  time  is  determined  by 
an  interpolation  procedure  in  the  program.  The  components  of  accelera- 
tion can  be  expressed  in  either  a launch-centered  inertial  or  launch- 
centered  rotating  coordinate  system,  depending  on  which  system  the 
inertial  platform  uses. 

When  the  missile  is  in  free  fall  under  the  Influence  of  gravity, 
acceleration  is  determined  analytically,  with  both  drag  and  thrust  set 
equal  to  zero. 

Earth  Model 

The  earth  can  be  chosen  to  be  either  rotating  or  nonrotating  and 
either  oblate  or  spherical.  If  an  oblate  earth  is  chosen,  the  term 

is  included  in  the  potential  function  and  the  difference  between  geodetic 
and  geocentric  latitude  is  taken  into  account. 
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TRACKER  CHARACTERISTICS 
Number  o f_  Tr acker s 

In  its  present  form,  the  program  can  handle  a maximum  of  24 
trackers,  each  of  which  can  make  a maximum  of  four  different  types  of 
measurement  at  a given  observation  time. 

Tracker  Location 

Each  of  the  24  trackers  can  be  located  on  the  ground,  on  an  air- 
craft, or  on  a satellite.  For  a ground-based  tracker,  it  is  necessary 
to  specify  its  latitude,  longitude,  and  altitude  relative  to  the  ro- 
tating earth.  It  is  assumed  that  the  aircraft  carrying  trackers  are 
flying  great  circle  routes  relative  to  the  rotating  earth.  Thus,  it 
is  necessary  to  specify  not  only  their  latitude,  longitude,  and  alti- 
tude, but  also  their  velocity  and  heading  as  well  as  the  time  at  which 
these  conditions  exist.  From  these  initial  position  data,  the  position 
and  velocity  at  some  future  time  can  be  determined.  The  satellites 
bearing  trackers  are  on  either  circular  or  elliptical  orbits  fixed  in 
inertial  space.  As  in  the  case  of  aircraft,  it  is  necessary  to  specify 
their  latitude,  longitude,  and  altitude  relative  to  the  rotating  earth 
at  a given  time,  as  well  as  the  velocity,  path  elevation  angle,  and 
azimuth  angle  relative  to  inertial  space.  From  these  initial  data, 
the  position  and  velocity  of  the  satellite  at  any  other  time  can  be 
determined . 

A fourth  type  of  tracker  which  can  be  introduced  is  an  altimeter 
located  in  the  missile.  The  most  convenient  way  to  implement  this  in 
the  program  is  as  a hypothetical  shadow  tracker  that  moves  along  the 
surface  of  the  earth  below  the  missile  and  measures  the  range  to  the 
missile.  These  range  measurements  are  identical  with  the  altitude 
measurements  that  would  be  produced  by  an  altimeter  in  the  missile. 

It  should  be  noted  that  these  altitude  measurements  are  assumed  to  be 
relative  to  an  ellipsoidal  earth,  since  the  program  does  not  model 
local  altitude  variations  of  the  earth's  surface. 

It  is  also  possible  to  introduce  biases  in  the  latitude,  longitude, 
altitude,  and  velocity  of  the  initial  positions  of  the  trackers — with 
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the  exception  of  the  shadow  trackers.  These  biases  can  be  specified 
directly  or  can  be  selected  from  gaussian  distributions  whose  standard 
deviations  have  been  specified. 

Tracker  Limitations 

There  are  a number  of  restrictions  to  the  operation  of  the  trackers 
which  can  be  set  individually  for  each  of  the  24  trackers.  The  first 
of  these  is  the  time  interval  during  which  the  tracker  is  in  operation 
and  available  to  make  measurements.  In  the  program,  it  is  possible  to 
specify  a maximum  of  four  time  intervals  during  which  a given  tracker 
is  operational. 

It  Is  also  possible  to  set  the  minimum  interval  between  observa- 
tions for  each  tracker.  It  should  be  noted  that,  if  a given  tracker 
is  capable  of  more  than  one  type  of  measurement,  all  of  these  measure- 
ments are  assumed  to  be  made  simultaneously  and  the  interval  specified 
above  represents  the  minimum  time  before  another  set  of  such  measurements 
can  be  made  by  the  same  tracker. 

Finally,  the  ability  of  a tracker  to  observe  a target  may  be  re- 
stricted by  the  target  elevation  angle  relative  to  the  tracker.  This 
can  be  specified  in  the  form  of  a minimum  and  maximum  permissible  ele- 
vation angle  for  each  tracker. 

MEASUREMENT  CHARACTERISTICS 

Types  of  Measurement 

There  is  a maximum  of  six  different  types  of  measurements  which 
can  be  made  by  any  one  tracker:  range,  range  rate,  azimuth  angle, 

elevation  angle,  hyperbolic  range,  and  hyperbolic  rate.  The  details 
of  the  various  types  of  measurement  are  described  in  Sec.  IV.  Although 
six  types  are  available,  it  is  not  anticipated  that  more  than  four 
would  be  used  by  any  one  tracker  so  that,  with  24  trackers,  the  maximum 
number  of  measurements  at  any  one  time  would  be  96. 

Measurement  Biases 


A bias  can  be  specified  for  each  type  of  measurement  and  each 
tracker  that  need  not  be  the  same  for  all  trackers.  The  magnitudes  of 


these  biases  can  be  inpat  directly  or  selected  from  gaussian  distri- 
butions with  appropriate  standard  deviations. 

Standard  Deviations 

For  each  type  of  measurement  which  can  be  made  by  a given  tracker, 
it  is  necessary  to  assign  a nonzero  standard  deviation  that  can  vary 
from  tracker  to  tracker. 
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IV.  MEASUREMENT  PROCEDURE 


This  section  describes  the  mechanics  of  the  measurement  procedure 
at  each  measurement  time. 


TRACKER  SELECTION 

Each  of  the  trackers  is  examined  in  sequence  to  determine  whether 
it  is  available  to  make  an  observation  at  a given  measurement  time. 

If  this  condition  is  satisfied,  the  position  and  velocity  of  the  nomi- 
nal missile  relative  to  this  tracker  are  determined  in  tracker  coordi- 
nates, represented  by  x^,  yR,  zR,  yR,  and  zR.  Prom  these  coordinates, 
the  missile’s  range,  range  rate,  azimuth  angle,  and  elevation  angle  rela- 
tive to  the  tracker  are  determined  from  the  following  relations: 


Range 


R 


(26) 


Range  rate 


VR  + yRyR 


ZRZR 


(27) 


Elevation  angle 


Azimuth  angle 


(28) 


(29) 


At  this  point,  a check  is  made  to  be  sure  that  the  elevation  angle  Y 
lies  between  the  maximum  and  minimum  values  for  the  tracker.  In  addi- 
tion, if  the  tracker  is  on  a satellite,  a check  is  made  to  determine 
whether  the  line  of  sight  between  the  tracker  and  missile  is  unob- 
structed. If  either  of  these  conditions  is  violated,  the  tracker  is 
not  available  and  a new  one,  if  available,  must  be  selected. 


-24- 


After  the  above  visibility  tests  have  been  satisfied,  the  biased 

values  of  the  variables  determined  in  Eqs.  (26)  through  (29)  are  found 

by  computing  the  position  and  velocity  of  the  nominal  missile  relative 

to  the  biased  tracker  position  in  tracker  coordinates  represented  by 

_ v , i , i , and  znt.  The  biased  versions  of  Eqs.  (26) 

*RB’  yRB*  ZRB  *RB*  yRB’  R6 

through  (29)  have  the  form 


SRU 


Dl 

Biased  range  ^ " *B0  + 6RB  + 0.1  + 0. 9 sin 


(30) 


Biased 
range  rate 


• *RB*RB  * yRByRB 

rbo 


+ ZRBZRB 


+ 


+ 


0.1  + 0.9  sin  Yjjq 


(31) 


Biased  ^yBl 

elevation  YB  - YB0  + 6YB  + 0. 1 + 0.9  sin  Y0O 

angle 


Biased 

azimuth 

angle 


(32) 


(33) 


where 


R 


B0 


and 


rB0 


tan 


-1 


RB 


4 


2,2 

*RB  + yRB- 


The  quantities  6Rb>  6Rb,  6B0,  and  <5y0  are  the  fixed  measurement  biases 
for  the  given  tracker,  whereas  the  terms  involving  6R0r  «Rb1»  and  6yBl 
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represent  the  measurement  biases  due  to  tropospheric  effect*!  that  are 
elevation  angle  dependent. 

MEASUREMENT  TYPES 

The  computations  indicated  thus  far  in  this  section  are  made  for 
each  available  tracker,  regardless  of  the  type  of  measurements  it  is 
capable  of  making.  In  the  succeeding  paragraphs  describing  the  various 
measurement  procedures,  the  computations  are  made  only  if  the  tracker 
is  capable  of  the  specified  type  of  measurement. 

Range  Measurement 

If  the  tracker  is  intended  to  measure  range,  the  variance  of  the 
measurement  is  determined  from  the  relation 

22  ^7  ^19 

°R  " C1  + 2 + CHR  » (34) 

R 1 sin  y 1J 

where  the  first  term  is  due  to  random  time  delays  in  the  measuring 
system,  the  second  to  tropospheric  refraction,  and  the  third  to  the 
slgnal-to-nolse  ratio.  Since  the  model  used  to  compute  the  tropo- 
spheric term  assumes  that  the  tracker  is  on  the  ground  looking  up  at 
the  missile,  it  is  probably  desirable  to  omit  this  term  when  dealing 
with  satellite-borne  trackers  by  setting  equal  to  zero. 

The  reference  range  R^  between  the  nominal  tracker  and  missile 
positions  is  now  set  equal  to  the  value  of  R computed  in  Eq.  (26),  and 
the  measured  range  is  given  by 

Rjj  - Rb  + 6R  , (35) 

where  R^  is  obtained  from  Eq,  (30)  and  dR  is  a random  variation  in 
range  selected  from  a gauss ian  distribution  with  a standard  deviation 
aR  given  by  Eq.  (34). 


I 


f ! 


i 

i i 


t 

l 

V . 

p ’ 

p- 

f, 

% ' 
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The  H matrix  row  associated  with  this  range  measurement  is  com- 
puted by  first  finding  the  B matrix  row,  which  is  given  as 


0 


, , ■ 3R  SR  SR  SR  SR  SR 

11  »2r3*r3*r3*rJ 


>*R  *yR 

whose  elements  can  be  evaluated  by  means  of  Eq.  (26)  as 


Sr  \ 


SR  y_t 
3yR‘R 


88  3 

82r-8 


88  „ 
"k  “ 

88  r. 

-T—  « 0 

SyR 

Si  u 


The  H matrix  row  is  then  determined  by  the  relation 

[H]  * {B]  x [A]  , 


(36) 


(37) 


(38) 


where  the  A matrix,  which  is  derived  in  the  appendix,  has  the  form 


3xr 

9xr 

9XR 

9xr 

9*R 

9*R 

3XX 

ax2 

ax3 

3X4 

ax. 

D 

3X6 

3yR 

3yR 

3yR 

s*r 

3yR 

3yR 

3XX 

3X, 

4_ 

ax3 

ax4 

ax5 

9X6 

3zr 

3zr 

3zr 

3 ZR 

3zr 

3 ZR 

ax3 

3X2 

3X3 

ax4 

ax5 

ax6 

3*r 

3*R 

3*R 

3*r 

3*r 

3*r 

3xi 

3X2 

ax3 

ax4 

3X 

SX6 

3yR 

3yR 

3yR 

3yR 

9*R 

3yR 

3X1 

3X2 

3X3 

ax4 

ax5 

3Xb 

3zR 

3zr 

3zr 

3zr 

3zr 

3zr 

3X1 

ax., 

4. 

ax3 

ax. 

4 

3X5 

3X6 

Thus,  Eq.  (38)  gives  the  H matrix  in  the  form  shown  in  Eq.  (22),  where 
Z is  replaced  by  R as  the  measured  quantity. 


Range  Rate  Measurement 

If  the  tracker  is  intended  to  measure  range  rate,  the  variance  of 
the  measurement  is  determined  from  the  relation 


O?  - C2  + f—  + C ,R  20  (40) 

R 1 sin  y 


in  a form  similar  to  that  used  for  the  range  variance  in  Eq.  (34). 
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The  reference  range  rate  Rq  between  the  nominal  tracker  and  mis- 
sile positions  is  set  equal  to  R,  computed  in  Eq.  (27),  and  the  mea- 
sured range  rate  R,.  is  given  by 


«M  - *8  + • 


(41) 


• * 
where  Rg  is  obtained  from  Eq.  (31)  and  6R  is  a random  variation  in 

range  rate  which  can  be  selected  from  a gaussian  distribution  with  a 

standard  deviation  given  by  Eq.  (40). 

As  in  the  case  of  range  measurement,  another  row  of  the  B matrix 

is  determined  in  the  form 


[B] 


3R  3_R_  3R__  3R_  3R_  3R_ 

3X8  % 3zR  3iR  3*R  3zR  ’ 

L -J 


(42) 


where  the  elements  of  this  row  matrix  can  be  determined  from  Eq.  (27) 
in  the  form 

3R  _ *R  ^R 

3xr  = r "r2 


9yR  R R2 


3R  _ *R 
3Xr  = r 


3R 

9yR"R 


3R  ^R 
3zR*R  * 


(43) 
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The  corresponding  row  of  the  H matrix  is  again  obtained  by  means  of 
Eq.  (38). 

Elevation  Angle  Measurement 

If  the  tracker  is  intended  to  measure  elevation  angle,  the  vari- 
ance of  the  measurement  is  determined  from  the  relation 


, _ cos  Y . q p 

y C3  + C9  . 2 + C15R 

' sin  y 


where  the  tropospheric  term  has  a somewhat  different  dependence  on 
the  elevation  angle  than  in  the  case  of  the  range  and  range  rate 
measurements . 

The  reference  elevation  angle  Yc  of  the  nominal  missile  relative 
to  the  nominal  tracker  is  set  equal  to  the  value  of  y computed  in  Eq. 
(28).  The  measured  elevation  angle  YM  is  given  by 


yM  ° + 4Y  • 


where  Yn  is  given  by  Eq.  (32)  and  6y  is  a random  elevation  angle  vari- 
B 

ation  selected  from  a gaussian  distribution  with  a standard  deviation 
0 given  by  Eq.  (44). 

The  row  of  the  B matrix  corresponding  to  this  measurement  has  the 


|"lX_  lL_  2L_  13L.  1 

K 3zR  8iR  »R  3iR  J ’ 


where  the  elements  of  this  row  matrix  can  be  determined  from  Eq.  (28) 
in  the  form 


The  corresponding  row  of  the  H matrix  is  again  obtained  by  means 
of  Eq.  (38). 


Azimuth  Angle  Measurement 

If  the  tracker  is  intended  to  measure  azimuth  angle,  the  variance 
of  this  measurement  is  given  by 


2 

aB 


2 

C,  + C1tR 
A 16 


(48) 
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The  reference  azimuth  angle  6 of  the  nominal  missile  relative 

v 

to  the  nominal  tracker  is  set  equal  to  the  value  of  6 obtained  from 
Eq.  (29);  the  measured  azimuth  angle  8w  is  given  by 

eM  - SB  + 68  , (49) 

where  8„  is  given  by  Eq.  (33)  and  63  is  a random  azimuth  angle  varia- 

D 

tion  which  can  be  selected  from  a gaussian  distribution  having  a stand- 
ard deviation  given  by  Eq.  (48). 

The  row  of  the  B matrix  corresponding  to  this  measurement  has  the 

form 


r fas  33  38  98  as  as 

[ J = L3xR  *R  3zR  3*R  3*R  3zR  J ’ 

where  the  elements  of  this  row  can  be  determined  from  Eq.  (29)  as 
follows: 


93  _ *R 


98 

3yR 


9J_ 

3zr 


0 


98 

3*R 


0 


(51) 
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The  corresponding  row  of  the  H matrix  is  determined  by  means  of 
Eq.  (38). 

Hyperbolic  Range  Measurement 

This  measuring  system  involves  a one-way  range  measurement  in 
which  each  tracker  radiates  a continuous  signal  with  a modulation  that 
Identifies  its  time  of  origin  when  received  by  the  missile.  At  a given 
measurement  time,  the  missile  can  determine  the  time  of  propagation  from 
each  of  the  observed  trackers  from  the  difference  between  the  measure- 
ment time  as  measured  by  the  onboard  clock  and  the  respective  times  of 
origin  as  determined  by  a ground-based  clock.  If  both  of  these  clocks 
are  synchronized,  the  resulting  propagation  times  can  be  ised  to  de- 
termine the  instantaneous  ranges  to  each  of  the  observing  trackers. 

Since  there  is  usually  a drift  between  the  two  clocks,  it  is  preferable 
to  select  one  of  the  trackers  as  a reference  and  determine  the  differ- 
ences in  propagation  time  between  each  of  the  other  trackers  and  the 
selected  reference  tracker.  In  this  way,  the  error  in  the  onboard 
clock  can  be  eliminated. 

The  implementation  of  this  method  is  Initially  identical  with  that 
described  for  the  two-way  range  measurement.  The  standard  deviation  of 
the  measurement  is  given  by 

°HR  - C5  + - V + ex/”  • <52> 

HR  5 Bln  Y 

which  is  of  the  same  form  as  Eq.  (34)  but  differs  in  the  values  of  C's 
because  of  the  one-way  nature  of  the  measurement.  The  reference  range 
Rc,  the  measured  range  R^,  the  elements  of  the  B matrix,  and  the  cor- 
responding elements  of  the  H matrix  are  determined  for  each  tracker  as 
before. 

To  use  the  hyperbolic  range  method,  at  least  two  trackers  must  be 
observed  at  a given  measurement  time.  If  this  condition  is  satisfied, 
the  tracker  with  the  minimum  range  is  selected  as  a reference.  A range 
difference  measurement  can  then  be  determined  for  each  of  the  other 
trackers  from  the  difference  of  its  measured  range  and  that  of  the 
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reference  tracker.  Similarly,  the  nominal  range  difference  for  each 
tracker  is  its  nominal  range  minus  that  of  the  reference  tracker. 
Finally,  the  elements  of  the  H matrix  for  these  range  difference  mea- 
surements are  obtained  from  the  difference  of  the  corresponding  ele- 
ments of  the  H matrix  computed  above  for  the  given  tracker  and  those 
for  the  reference  tracker  in  the  form 


(53) 


Thus,  the  number  of  independent  hyperbolic  range  measurements  is  one 

less  than  the  number  of  observing  trackers. 

Another  version  of  the  hyperbolic  range  method  eliminates  the 

reduction  in  the  number  of  measurements  at  each  observation  time  after 

the  first  by  saving  the  minimum  range  measurement  and  the  nominal  range 

from  the  previous  observation  time,  as  well  as  the  corresponding  row 

of  the  H matrix,  H. (R  . ).  This  H matrix  row  is  propagated  forward  to 
1 min 

the  new  observation  time  by  the  relation 

H(Rref>  ' H(W  l+l'1  • 


where  [<j>)  is  the  upper  left  hand  6 x 6 of  the  state  transition  matrix 
$.  The  saved  values  of  the  range  measurement,  the  nominal  range,  and 
the  propagated  H matrix  elements  can  now  be  used  as  the  reference  values 
tc  be  subtracted  out  at  the  new  observing  time.  Since  these  values  can 
also  be  subtracted  from  the  observations  of  the  tracker  which  provided 
the  previous  minimum  range,  the  number  of  independent  hyperbolic  range 
measurements  is  now  equal  to  the  number  of  observing  trackers. 


-34- 


A final  version  of  the  hyperbolic  range  measurement  method  uses 
a different  reference  for  each  tracker — the  previous  measurement  made 
by  the  given  tracker.  All  of  the  range  measurements  and  nominal 
ranges  from  the  previous  observation  time  must  thus  be  saved,  as  well 
as  the  entire  H matrix,  which  is  then  propagated  forward  to  the  present 
time  by  means  of  Eq.  (54).  The  resulting  matrix,  together  with  the 
saved  measured  ranges  and  nominal  ranges,  can  be  subtracted  from  the 
corresponding  data  for  the  respective  trackers  to  provide  the  desired 
hyperbolic  range  data.  It  is  seen  that,  in  this  case,  a hyperbolic 
range  measurement  can  only  be  made  by  a tracker  if  it  was  observed  at 
the  previous  observation  time  as  well  as  at  present.  This  approach 
has  the  advantage  that  it  tends  to  cancel  out  the  individual  tracker 
position  biases  as  well  as  the  clock  error;  it  is  referred  to  as  the 
delta  range  method. 

Hyperbolic  Range  Rate  Measurement 

This  measuring  system  is  entirely  analogous  to  the  hyperbolic 
range  method  described  above;  it  has  a standard  deviation  of  the  form 
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In  this  case,  the  range  rate  for  each  tracker  is  determined  onboard 
the  missile  by  comparing  the  received  frequency  with  the  frequency  of 
an  onboard  oscillator  which  presumably  is  matched  to  the  radiated 
frequency.  To  eliminate  the  frequency  error  between  the  onboard  os- 
cillator and  the  radiated  frequency,  a differencing  technique  is  used 
in  which  the  data  corresponding  to  the  minimum  range  tracker  are  sub- 
tracted from  the  corresponding  data  for  the  other  trackers  to  give  the 
desired  hyperbolic  range  rate  measurements  in  a form  appropriate  for 
processing  by  the  Kalman  filter.  This  method  also  has  the  option  of 
establishing  a reference  set  of  data  by  propagating  the  data  from  the 
minimum  range  tracker  at  the  previous  measuring  time  to  the  present 
time  by  means  of  the  state  transition  matrix,  in  accordance  with 


Eq.  (54).  Although  it  could  be  done,  the  program  does  not  provide  for 
saving  individual  measurements  and  H matrix  elements  for  all  trackers 
as  references  for  the  next  measurements  by  these  systems,  as  was  done 
in  the  delta  range  case. 
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V.  RTRAJ  PROGRAM  OPERATION 


INPUT  DECK 

With  the  exception  of  one  card,  the  RTRAJ  input  deck  uses  the 
FORTRAN  namelist  method  of  data  input.  This  subsection  lists  the 
input  variables  with  their  definitions  and  their  default  values  in 
the  event  they  are  not  specified.  The  namelists  and  the  one  formatted 
card  input  appear  in  their  required  order,  although  the  order  of  the 
variables  within  a given  namelist  is  immaterial. 

Namelist /BIASIN/ 

IRUN  Number  of  missile  flights  submitted  in  this  input 

deck;  a default  value  of  one. 

SPRANG  Standard  deviation  in  range,  in  feet,  of  the  gaussian 
distribution  from  which  the  range  measurement  biases 
of  the  various  trackers  can  be  selected. 

SDRRT  Standard  deviation  in  range  rate,  in  feet  per  second, 

of  the  gaussian  distribution  from  which  the  range  rate 
measurement  biases  of  the  various  trackers  can  be 
selected, 

SDELEV  Standard  deviation  in  elevation  angle,  in  milliradians , 
of  the  gaussian  distribution  from  which  the  elevation 
angle  measurement  biases  of  the  various  trackers  can 
be  selected. 

SDAZIM  Standard  deviation  in  azimuth  angle,  in  milliradians, 
of  the  gaussian  distribution  from  which  the  azimuth 
angle  measurement  biases  of  the  various  trackers  can 
be  selected. 

SDPNG  Standard  deviation  in  range,  in  feet,  of  the  gaussian 

distribution  from  which  the  tropospheric  range  measure- 
ment bias  at  90-deg  elevation  angle  for  the  various 
trackers  can  be  selected, 

SPRAT  Standard  deviation  in  range  rate,  In  feet  per  second, 

of  the  gaussian  distribution  from  which  the  tropospheric 
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range  rate  measurement  bias  at  90-deg  elevation  angle 
for  the  various  trackers  can  be  selected. 

SDJE1.V  Standard  deviation  in  elevation  angle,  in  milliradians, 

of  the  gaussian  distribution  from  which  the  tropo- 
spheric elevation  angle  measurement  bias  at  90-deg 
elevation  angle  for  the  various  trackers  can  be  selected. 

SPLAT  Standard  deviation  in  latitude,  in  feet,  of  the  gauss "an 

distribution  from  which  the  latitude  biases  of  the  vari- 
ous trackers  can  be  selected. 

SDLON  Standard  deviation  in  longitude,  in  feet,  of  the  gaussian 

distribution  from  which  the  longitude  biases  of  the  vari- 
ous trackers  can  be  selected. 

SDALT  Standard  deviation  in  altitude,  in  feet,  of  the  gaussian 

distribution  from  which  the  altitude  biases  of  the  vari- 
ous trackers  can  be  selected. 

SDVEL  Standard  deviation  in  tracker  velocity,  in  feet  per 

second,  of  the  gaussian  distribution  from  which  the 
biases  in  velocity  of  the  various  trackers  can  be  selected. 

SX(I)  The  six  standard  deviations  in  missile  position,  in  feet 

(I  * 1,2,3),  and  velocity,  in  feet  per  second  (1  = 4,5,6), 
of  the  gaussian  distributions  from  which  missile  position 
and  velocity  biases  can  be  selected. 

The  default  values  of  all  these  standard  deviations  are  zero. 

Formatted  Card  Input 

The  format  for  this  lone  input  card  is  D10.0,  15,  15,  and  the 
variables  ure  as  follows: 

PRNTDT  Interval  in  flight  time,  in  seconds,  at  which  output 
data  are  to  be  printed. 

ISTG  A flag  to  indicate  whether  this  stage  is  the  last  one 

of  this  flight. 

0 This  is  the  final  stage  of  this  flight. 

1 Another  stage  follows. 
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IBST  A flag  to  indicate  the  type  of  acceleration  profile 

used;  it  can  have  an  integral  value  from  0 to  2 with 
the  following  meanings: 

0 The  missile  is  in  free  fall  and  its  accelera- 
tion is  determined  analytically. 

1 The  missile  acceleration  is  input  in  tabular 
form  as  a function  of  time  in  launch-centered 
inertial  coordinates.  It  also  implies  that 
the  inertial  platform  associated  with  the  1MU 
is  stabilized  in  this  same  coordinate  system. 

2 The  missile  acceleration  is  input  in  tabular 
form  as  a function  of  time  in  launch-centered 
rotating  coordinates.  It  also  implies  that 
the  inertial  platform  is  stabilized  in  this 
same  coordinate  system. 


Name  list /STATIN/ 


IN?Q 


INPX 


A flag  to  describe  the  form  of  the  input  of  the  ini- 
tial covariance  matrix  Q;  it  can  have  integral  values 
between  -4  and  +4.  If  the  sign  of  INPQ  is  negative, 
only  the  diagonal  elements  of  the  covariance  matrix  are 
entered,  as  DATQ;  a positive  sign  means  that  the  entire 
covariance  matrix  is  entered,  as  Q.  The  magnitude  of 
INPQ  indicates  the  coordinate  system  in  which  the  co- 
variance  matrix  is  expressed. 

0,1  Earth-centered  inertial. 

2 In-plane  inertial. 

3 Earth-centered  rotating. 

4 In-plane  rotating. 

INPQ  should  only  be  specified  in  the  first  stage  of 
the  flight  and  then  only  if  it  differs  from  the  default 
value  of  zero. 

A flag  to  determine  the  coordinate  system  in  which  the 
estimated  state  vector  DX  is  expressed.  It  can  have  an 
integral  value  from  0 to  4 with  the  same  meanings  as 
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IX 


m 


IS 


those  for  the  magnitude  of  INPQ.  As  above,  it  is  only 
necessary  to  specify  INPX  in  the  first  stage  of  the 
flight  and  then  only  if  it  differs  from  the  default 
value  of  zero. 

This  flag  controls  the  printout  of  the  state  variables 
X and  the  state  vector  DX  and  can  have  the  integral 
values  of  0 or  1 with  the  meanings: 

0 Do  not  print  X and  DX. 

1 Print  X and  DX. 

The  default  value  of  IX  is  zero. 

A flag  to  control  the  printout  of  the  covariance  matrix 
as  well  as  the  associated  eigen  values  and  eigen  vectors. 
It  can  have  integral  values  from  0 to  7 with  the  follow- 
ing meanings: 

0,2, 4, 6 Do  not  print  the  upper  left  6 x 6 of  the 
covariance  matrix. 

0,1, 4, 5 Do  not  print  the  remaining  tridiagonal 
of  the  covariance  matrix  which  includes 
those  elements  in  rows  7 through  51  from 
column  one,  out  to  and  including  the 
diagonal . 

0,1, 2, 3 Do  not  print  the  eigen  values  and  eigen 
vectors  associated  with  the  covariance 
matrix. 

7 Print  everything. 

It  is  seen  that,  by  a suitable  choice  of  IQ,  any  com- 
bination of  the  three  items  listed  above  can  be  printed 
out.  The  default  value  of  IQ  is  zero. 

A flag  to  determine  the  coordinate  system  in  which  the 
state  vector,  the  covariance  matrix,  and  the  eigen  vec- 
tors will  be  printed  out.  It  can  have  integral  values 
from  0 through  4 with  the  following  meanings: 

■ 1 Printout  is  in  earth-centered  inertial 
coordinates. 

2 Printout  is  in  in-plane  inertial  coordinates. 
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DX(I) 


DATQ(I) 


Q(l.J) 


IMODE 


3 Printout  is  in  earth-centered  rotating 
coordinates . 

4 Printout  is  in  in-plane  rotating  coordinates. 
The  default  value  of  IS  is  zero. 

A 51-element  array  comprising  the  initial  values  of 
the  estimated  state  vector,  where  the  first  six  asso- 
ciated with  the  missile’s  position  and  velocity  are 
expressed  in  the  coordinate  system  specified  by  INPX, 
and  the  remaining  45  associated  with  the  IMU  are  inde- 
pendent of  the  coordinate  system.  The  default  values 
for  all  of  the  DX's  are  zero. 

A 51-element  array  specifying  the  initial  diagonal 
elements  of  the  covariance  matrix  that  are  the  variances 
of  the  state  variables  expressed  in  the  coordim  te  sys- 
tem specified  by  INPQ.  This  array  is  input  only  when 
the  value  of  INPQ  is  less  than  zero.  The  default  value 
of  DATQ(l)  is  zero. 

A 51  * 51  array  containing  the  initial  elements  of  the 
covariance  matrix  that  are  entered  only  when  INPQ  is 
greater  than  zero.  If  Q(I,J)  is  entered  when  INPQ  is 
less  than  zero,  its  diagonal  elements  will  be  super- 
seded by  the  values  of  DATQ.  The  default  values  for 
all  of  the  elements  of  Q(1,J)  are  zero. 

A flag  to  determine  which  mode  of  operation  of  the  pro- 
gram is  in  use,  as  described  in  Sec.  Ill  on  Program 
Options.  It  can  have  the  values  0,  1,  or  2 with  the 
following  meanings: 

0 The  program  operates  in  Mode  0,  in  which  co- 
variance  is  propagated  along  the  trajectory 
and  updated  by  the  Kalman  filter. 

1 The  program  operates  in  Mode  1,  in  which  ob- 
servational data  are  generated  and  recorded 
for  later  use  in  Mode  2. 

2 The  program  operates  in  Mode  2,  in  which  a de- 
termination of  the  state  vector  is  made  using 
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observational  data  generated  either  by  Mode  1 
or  from  the  tracking  of  an  actual  missile. 


Namelist /TRNFIN/ 

IOBS  A flag  to  control  the  printout  of  the  tracker  obser- 
vational data,  which  include:  the  tracker  identifi- 

cation numbers;  its  latitude,  longitude,  and  altitude; 
its  path  elevation  angle,  azimuth  angle,  and  speed; 
and,  in  addition,  the  magnitude  of  the  measured  quan- 
tity and  the  number  of  that  type  of  measurement  that 
has  been  made  by  that  tracker.  This  flag  is  an  inte- 
ger with  the  following  interpretations: 

0 Do  not  print  observational  information. 

>0  Print  observational  data  only  at  the  interval 
specified  by  PRNTDT.  Thus,  if  the  observa- 
tion interval  is  less  than  the  print  interval, 
only  those  observations  at  the  specified 
printout  time  are  printed. 

<0  Observational  information  for  all  measure- 
ments is  printed  out  regardless  of  the  print 
interval. 

The  default  value  of  I0BS  is  zero. 

IR  A flag  to  control  the  printout  of  the  standard  devia- 

tions of  the  measurements.  It  has  an  integer  value 
which  is  interpreted  as  follows: 

0 Do  not  print  standard  deviations  of  measure- 
ments . 

^0  Print  standard  deviations  of  measurements. 

The  default  value  of  IR  is  zero. 

IH  A flag  to  control  the  printout  of  the  H matrix  gen- 

erated by  the  measurements.  It  is  an  integer  which 
is  interpreted  as  follows: 

0 Do  not  print  the  H matrix. 

Print  the  H matrix. 

The  default  value  of  IH  is  zero. 
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TIN 


IRND 


GHA 


TMAX 


ICOMP 


DELMN 


The  starting  time  of  the  run,  in  seconds;  it  is  only 
entered  in  the  first  stage  since  it  is  automatically 
set  to  its  appropriate  value  in  succeeding  stages. 

The  default  value  of  TIN  is  zero. 

A flag  to  determine  whether  random  variations  are 
introduced  into  the  tracker  measurements.  If  IRND  is 
zero,  no  such  variations  are  included;  if  it  is  not 
zero,  then  a random  variation  is  selected  from  a gauss- 
ian  distribution  whose  standard  deviation  is  that 
specified  for  the  given  type  of  measurement.  This 
variation  is  then  added  to  the  already  biased  value 
of  the  measurement.  The  default  value  of  IRND  is  zero. 
The  initial  value  of  the  Greenwich  Hour  Angle,  which 
is  the  longitude  difference  between  Greenwich  and  the 
first  point  of  Aries  and  is  expressed  in  degrees.  It 
is  specified  only  for  the  first  stage,  since  it  is 
automatically  reset  to  its  appropriate  value  at  the 
beginning  of  each  succeeding  stage.  The  default  value 
of  GHA  is  zero. 

The  termination  time  of  this  stage  of  the  flight,  mea- 
sured in  seconds  from  the  start  of  the  flight.  If  this 
stage  is  intended  to  terminate  at  impact,  the  value  of 
TMAX  should  be  made  larger  than  the  estimated  time  to 
insure  that  the  missile  actually  reaches  impact.  The 
default  value  of  TMAX  is  zero. 

A flag  to  determine  whether  gravity  compensation  is 
used  in  the  IMU,  as  described  in  the  sections  on  gyro- 
scope and  accelerometer  errors.  If  ICOMP  is  a nonzero 
value,  compensation  is  used.  The  default  value  of 
ICOMP  is  zero. 

The  minimum  elevation  angle  of  a satellite  tracker 
relative  to  the  user  vehicle  such  that  the  tracker  can 
be  seen  by  the  user.  The  default  value  of  DELMN  is 


zero. 
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HMIN  If  a minimum  altitude  exists  along  the  path  connect- 

ing a satellite  tracker  and  a user  vehicle,  this  mini- 
mum cannot  be  less  than  HMIN  if  tracking  is  to  be 
accomplished.  HMIN  is  expressed  in  feet  and  has  a 
default  value  of  zero. 


Namelist/MODIN/ 

DTP  The  maximum  allowable  integration  interval,  in  seconds. 

IPRNT  A flag  to  control  the  printout  of  the  position  and 

velocity  of  the  missile.  If  its  value  is  less  than 
zero,  the  position  and  velocity  are  printed  out  at  the 
end  of  each  integration  interval;  if  IPRNT  is  greater 
than  zero,  these  data  are  printed  out  only  at  the 
specified  print  time.  The  absolute  value  of  IPRNT 
determines  the  coordinate  system  and  format  of  the 
position  and  velocity  printout  as  follows: 

0 Do  not  print  missile  position  and  velocity 
data. 

1 Print  missile  position  and  velocity  relative 
to  an  earth-centered  rotating  coordinate  sys- 
tem as  Xj^ , yR,  zR,  XR , yR,  and  zR,  measured  in 
feet  and  feet  per  second. 

2 Print  missile  position  and  velocity  relative 
to  an  in-plane  rotating  coordinate  system,  as 
latitude,  longitude,  altitude  (in  feet),  path 
elevation  angle  (in  degrees)  , path  azimuth 
angle  (in  degrees),  and  velocity  (in  feet  per 
second) . 

3 Print  missile  position  and  velocity  relative 
to  an  earth-centered  inertial  coordinate  sys- 
tem as  Xj,  y^,  Zj,  ij,  y^,  and  z^,  in  feet  and 
feet  per  second. 


The  default  value  of  IPRNT  is  zero. 
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A flag  to  determine  the  coordinate  system  in  which 
the  initial  position  and  velocity  of  the  missile  are 
expressed,  with  the  following  options: 

1 Earth-centered  rotating  coordinate  system. 

2 Rotating  topographical  latitude  and  longitude 
coordinate  system. 

3 Earth-centered  inertial  coordinate  system. 

4 Inertial  topographical  latitude  and  longitude 
coordinate  system. 

There  is  no  default  value  of  IP  and  it  is  entered 
only  for  the  first  stage. 

A three-elemer.t  array  that  specifies  the  initial  posi- 
tion of  the  missile;  however,  its  form  depends  on  the 
value  of  IP,  defined  above. 

• If  IP  equals  one,  the  three  elements  of  POS 
are  the  cartesian  coordinates  of  the  missile's 
position,  in  feet,  relative  to  an  earth- 
centered  rotating  coordinate  system. 

• If  IP  equals  two,  the  three  elements  of  POS 
are  the  latitude  (in  degrees) , the  longitude 
(in  degrees),  and  the  altitude  (in  feet)  of 
the  missile  relative  to  the  rotating  earth. 

• If  IP  equals  three,  the  three  elements  of  POS 
are  the  three  cartesian  coordinates,  in  feet, 
of  the  missile  relative  to  an  earth-centered 
inertial  coordinate  system. 

• If  IP  equals  four,  the  three  elements  of  POS 
are  the  latitude  (in  degrees),  the  longitude 
(in  degrees),  and  the  altitude  (in  feet)  of 
the  missile  relative  to  a nonrotating  earth. 

There  is  no  default  value  for  POS  and  it  should  only 
be  entered  in  the  first  stage  of  the  flight. 

A three-element  array  that  specifies  the  initial  ve- 
locity of  the  missile  although,  as  in  the  case  of  the 
POS  array,  its  form  depends  on  the  value  of  IP,  as 
follows : 
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TO 


TDOT(I) 


DDOT(I , 


• If  IP  equals  one,  the  three  elements  of  VEL 
are  the  cartesian  components  of  the  missile’s 
velocity,  in  feet  per  second,  relative  to  the 
earth-centered  rotating  coordinate  system. 

• If  IP  equals  two,  the  three  elements  of  VEL 
are  the  initial  path  elevation  angle  (in  de- 
grees), the  path  azimuth  angle  (in  degrees), 
and  the  velocity  (in  feet  per  second)  of  the 
missile  relative  to  the  rotating  earth. 

• If  IP  equals  three,  the  three  elements  of  VEL 
are  the  cartesian  components  of  the  initial 
missile  velocity,  in  feet  per  second,  relative 
to  the  earth-centered  inertial  coordinate 
system. 

• If  IP  equals  four,  the  three  elements  of  VEL 
are  the  path  elevation  angle  (in  degrees)  , 
the  path  azimuth  angle  (in  degrees),  and  the 
velocity  (in  feet  per  second)  relative  to  a 
nonrotating  earth. 

There  is  no  default  value  for  VEL  and  it  should  only 
be  entered  for  the  first  stage  of  the  flight. 

The  time  for  which  the  values  of  POS  and  VEL  are  de- 
fined; it  is  measured  in  seconds  relative  to  the  start 
of  the  run  and  has  no  default  value.  It  should  only 
be  specified  for  the  first  stage  of  the  flight. 

A 30-element  array  that  represents  the  values  of  time 
at  which  the  three  components  of  acceleration  are 
specified  when  a tabular  acceleration  input  is  used. 
There  is  no  default  value  for  TDOT  and  it  should  only 
be  entered  if  IBST  is  greater  than  zero.  The  final 
value  of  TDOT  should  be  the  same  as  TMAX  of  the  stage 
to  avoid  unspecified  values  of  input  acceleration. 

A 30  x 3 array  where  the  three  elements  of  the  1th 
row  are  the  three  components  of  missile  acceleration, 
in  feet  per  second  , at  time  TDOT(I).  These  components 
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PPSI 


IOBL 


IROT 


ALPHA(I) 


BETA ( I ) 


PHIA(I) 


are  in  launch-centered  inertial  coordinates  or  launch- 
centered  rotating  coordinates,  depending  on  whether 
1BST  is  equal  to  one  or  two.  There  is  no  default  value 
for  DDOT  and  it  should  only  be  entered  if  IBST  is 
greater  than  zero. 

The  launch  azimuth,  in  degrees,  measured  clockwise 
from  north.  It  has  no  default  value  and  must  be  spec- 
ified even  though  a path  azimuth  may  have  been  speci- 
fied as  VEL(2)  for  IP  aqual  to  two  or  four. 

A flag  to  specify  whether  the  earth  is  spherical  or 
oblate,  as  follows: 

0 Spherical  earth. 

1 Oblate  earth. 

The  default  value  is  zero. 

A flag  to  specify  whether  the  earth  is  nonrotating  or 
rotating,  as  follows: 

0 Nonrotatin?  earth. 

1 Rotating  earth. 

The  default  value  is  one. 

A 6-eleraent  array  representing  the  accelerometer  and 
gyroscope  orientation  angles  a^,  o^,  a^,  a^,  a^,  and 
a^,  as  defined  in  the  appendix.  The  default  values 
are  15,  -36,  -36,  -30,  -120,  and  -120  deg,  respectively. 
A 6-element  array  representing  the  accelerometer  and 
gyroscope  orientation  angles  B^,  ®3»  3^»  and 

8 y The  default  values  are  0,  28,  -28,  0,  0,  and 
90  deg,  respectively. 

A 6-element  array  representing  the  accelerometer  and 
gyroscope  orientation  angles  Y£>  Y3>  Y^>  Y^>  and 
yy  The  default  values  are  all  zero. 


Namelist/TKAKIN/ 

All  indexed  variables  in  this  namelist  have  a dimension  of  24, 
since  they  all  represent  characteristics  of  the  various  trackers  of 
which  there  can  be  a maximum  of  24. 
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ALA(I) 

ALO(I) 
ALT ( I ) 
S(I) 

V(I) 

B(I) 

PSI(I) 

GAM(I) 

TO  (I) 

TRTYP(I) 


ELMN(I) 


The  latitude  of  the  tracker,  in  degrees.  If  the  tracker 
Is  ground-based  or  air-based,  this  is  geodetic  latitude; 
otherwise,  it  is  geocentric  latitude.  The  default  value 
is  zero. 

The  longitude  of  the  tracker,  in  degrees.  The  default 
value  is  zero. 

The  altitude  of  the  tracker,  in  feet.  The  default  value 
is  zero. 

The  speed  of  an  airborne  tracker  relative  to  the  surface 
of  the  earth,  in  feet  per  second.  The  default  value  is 
zero. 

The  velocity  of  a satellite  tracker  relative  to  inertial 
space,  in  feet  per  second,  with  a default  value  of  zero. 
The  heading  of  an  airborne  tracker  relative  to  the  ro- 
tating earth,  in  degrees,  clockwise  from  north.  Its 
default  value  is  zero. 

The  heading  of  a satellite  tracker  in  inertial  space, 
in  degrees,  clockwise  from  north.  Its  default  value 
is  zero. 

The  path  elevation  angle  of  a satellite  tracker  in  in- 
ertial space,  in  degrees,  with  a default  value  of  zero. 
The  time  relative  to  the  start  of  the  run  at  which  the 
foregoing  variables  in  this  namelist  apply,  with  a de- 
fault value  of  zero. 

This  variable  is  alphanumeric  and  indicates  the  type 
of  tracker;  it  can  be  set  equal  to  the  following: 

• Ground. 

• Airplane. 

• Satellite. 

• Shadow. 

The  default  value  for  this  variable  is  BLANK. 

The  minimum  elevation  angle,  in  degrees,  of  the  target 
relative  to  the  tracker  for  satisfactory  observation. 

The  default  value  is  zero. 
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ELMX(I) 

ONA(I) 

OFFA(I) 

ONB(I) 

OFFB(I) 

PNC (I) 

OFFC(I) 

OND(l) 

OFFD(l) 

PELT (I) 

D1(I) 

D2(I) 

D3(I) 

DA(1) 

D5(I) 


The  maximum  elevation  angle,  in  degrees,  of  the  target 
relative  to  the  tracker  for  satisfactory  observation. 
The  default  value  is  100. 

The  starting  time  of  the  first  observing  interval  for 

the  tracker  measured,  in  seconds,  from  the  start  of 

20 

the  run.  The  default  value  is  -10 

The  terminal  time  of  the  first  observing  interval, 

20 

with  a default  value  of  +1 

The  starting  time  of  the  second  observing  interval, 

20 

with  a default  value  of  +10 

The  terminal  time  of  the  second  observing  interval, 

20 

with  a default  value  of  +10 

The  starting  time  of  the  third  observing  interval, 

20 

with  a default  value  of  +10 

The  terminal  time  of  the  third  observing  interval, 

20 

with  a default  value  of  +10 

The  starting  time  of  the  fourth  observing  interval, 

20 

with  a default  value  of  10 

The  terminal  time  of  the  fourth  observing  interval, 

20 

with  a default  value  of  10 

The  interval,  in  seconds,  between  observation  times 
for  the  tracker.  The  default  value  is  zero. 

The  constant  C^,  in  feet,  in  the  expression  for  the 
standard  deviation  in  the  range  measurement,  Eq.  (34). 
The  constant  C2,  in  feet  per  second,  in  the  expression 
for  the  standard  deviation  in  the  range  rate  measure- 
ment, Eq.  (40). 

The  constant  C^»  in  milliradians,  in  the  expression 
for  the  standard  deviation  of  the  elevation  angle 
measurement,  Eq.  (44). 

The  constant  C^,  in  milliradians,  in  the  expression 
for  the  standard  deviation  in  the  azimuth  angle  mea- 
surement, Eq.  (48). 

The  constant  C^,  in  feet,  in  the  expression  for  the 
standard  deviation  in  the  hyperbolic  range  measurement, 
Eq.  (52). 


— 
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D6(I) 


D7(I) 
D8(I) 
D9  (I) 


DRRT(I) 
DELEV (I) 

DAZIM(I) 

DRNG(I) 


The  constant  Cg,  in  feet  per  second,  in  the  expression 
for  the  standard  deviation  in  the  hyperbolic  range  rate 
expression,  Eq . (55). 

The  constant  in  Eq.  (34). 

The  constant  Cg  in  Eq.  (40). 

The  constant  Cg  in  Eq.  (44). 


Dll (I) 

The 

constant 

C11 

in 

Eq. 

(52). 

D12U0 

The 

constant 

C12 

in 

Eq. 

(55). 

D13(I) 

The 

constant 

C13 

in 

Eq. 

(34). 

D14 (I) 

The 

constant 

C14 

in 

Eq. 

(40). 

D15(I) 

The 

constant 

C15 

in 

Eq. 

(44). 

D16(I ) 

The 

constant 

C16 

in 

Eq. 

(48). 

D17  (I) 

The 

constant 

C17 

in 

Eq. 

(52). 

D18 (I ) 

The 

constant 

C18 

in 

Eq. 

(55). 

Dl 9 (I ) 

The 

constant 

C19 

in 

Eq. 

(34). 

D20(I) 

The 

constant 

C20 

in 

Eq. 

(40). 

D21(I) 

The 

constant 

C21 

in 

Eq. 

(44). 

D22  (I) 

The 

constant 

C22 

in 

Eq . 

(48). 

D23(I) 

The 

constant 

C23 

in 

Eq. 

(52). 

D24  (I) 
DRANG (I) 

The 

The 

constant  C„,  in  Eq. 
** 

bias  in  the  range 

(55). 

measurement 

The  bias  in  the  range  measurement,  in  feet  per  second. 
The  bias  in  the  elevation  angle  measurement,  in  milli- 
radians. 

The  bias  in  the  azimuth  angle  measurement,  in  milli- 
radians . 

The  tropospheric  range  measurement  bias,  in  feet,  for 
a 90-deg  elevation  angle. 


* 

The  default  values  for  the  constants  Dl  through  D24  are  zero; 
however,  it  is  imprrtant  that  the  standard  deviations  of  the  measure- 
ments as  given  by  Eqs.  (34),  (40),  (44),  (48),  (52),  and  (55)  should 
be  nonzero  for  each  type  of  measurement  made  by  a given  tracker.  This 
can  be  accomplished  by  assigning  nonzero  values  to  the  appropriate 
constants  in  these  expressions. 

** 

The  11  bias  terms  listed  above  have  default  values  of  zero,  and, 
if  these  biases  are  entered  in  the  namelist,  they  supersede  any  values 
that  mav  have  been  previously  selected  from  gaussian  distributions. 
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DRAT(I)  The  tropospheric  range  rate  measurement  bias,  in  feet 
per  second,  for  a 90-deg  elevation  angle. 

DELV(I)  The  tropospheric  elevation  angle  measurement  bias,  in 
milliradians,  for  a 90-deg  elevation  angle. 

DLAT(I)  The  bias  in  the  tracker  latitude,  in  feet. 

DLON(I)  The  bias  in  the  tracker  longitude,  in  feet. 

DALT(I)  The  bias  in  the  tracker  altitude,  in  feet. 

DVEL(I)  The  bias  in  the  tracker  velocity,  in  feet. 

NTR(I)  The  tracker  identification  number  associated  with  the 
tracker . 

AZIM(l)  If  set  equal  to  TRUE,  the  tracker  can  measure  azimuth 
* 

angle. 

ELEV (I)  If  set  equal  to  TRUE,  the  tracker  can  measure  elevation 
angle . 

RANGE (I ) If  set  equal  to  TRUE,  the  tracker  can  measure  range. 

RRT(I)  If  set  equal  to  TRUE,  the  tracker  can  measure  range 

rate. 

HRANGE(I)  If  set  equal  to  TRUE,  the  tracker  can  measure  hyper- 
bolic range. 

HRRT(I)  If  set  equal  to  TRUE,  the  tracker  can  measure  hyper- 
bolic range  rate. 

IPSMD  A flag  to  determine  the  part  cular  mode  of  hyperbolic 

range  measurement,  which  is  to  be  used  as  follows: 

0 The  minimum  range  tracker  is  used  as  the 
reference  at  each  observation  time. 

1 The  minimum  range  tracker  data  from  the  pre- 
vious measuring  time  are  propagated  forward 
using  the  state  transition  matrix,  to  be  used 
as  a reference  at  the  present  observation  time. 

2 The  previous  measurement  data  for  each  tracker 
are  propagated  forward  by  the  state  transition 
matrix,  to  be  used  as  references  for  the 


The  six  logical  variables  all  have  default  values  of  FALSE. 


Sti 
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respective  trackers  at  the  current  observation 
time  (delta  range  option). 

The  default  value  is  zero. 

IPRMD  A flag  to  determine  the  particular  mode  of  hyperbolic 

range  rate  measurement,  which  is  to  be  used  as  follows: 

0 The  minimum  range  tracker  is  used  as  a refer- 
ence at  each  observation  time. 

1 The  minimum  range  tracker  data  at  the  previous 
observation  time  are  propagated  forward,  using 
the  state  transition  matrix,  to  be  used  as  a 
reference  at  the  current  observation  time. 

Multistage  Operation 

The  organization  of  the  input  deck  described  above  is  that  re- 
quired for  the  first  stage  of  a missile's  flight.  If,  as  is  usually 
the  case,  there  is  more  than  one  stage  in  the  flight,  the  organization 
of  the  deck  is  as  follows: 

Stage  1 & BIASIN 

Unformatted  data  cards 

& END 

FORMATTED  DATA  CARD  (ISTG-1) 

& STATIN 

Unformatted  data  cards 

& END 
& TRNFIN 

Unformatted  data  cards 

& END 
& MODIN 

Unformatted  data  cards 

& END 
& TRAKIN 

Unformatted  data  cards 


& END 
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Stage  2 FORMATTED  DATA  CARD  (ISTG=0) 

& STATIN 

Unformatted  data  cards 

& END 
& TRNFIN 

Unformatted  data  cards 
& MOD IN 

Unformatted  data  cards 

& END 
& TRAKIN 

Unformatted  data  cards 

& END 

The  example  above  represents  only  two  stages,  but  the  same  oper- 
ation could  be  repeated  for  any  number  of  stages,  with  ISTG  set  equal 
to  unity  In  all  but  the  last  stage.  If  an  Input  variable  is  not 
specified  in  its  namelist,  it  will  retain  its  value  from  the  previous 
stage,  or  its  default  value  if  the  input  is  for  the  first  stage.  As 
is  frequently  the  case,  there  may  be  no  changes  in  the  values  of  the 
variables  in  a given  namelist  in  a succeeding  stage.  Nevertheless, 

It  is  still  necessary  to  include  the  namelist  title  card  and  END  card 
in  the  deck,  even  though  there  are  no  data  cards  between  them. 

Multiple  Runs 

It  is  also  possible  to  submit  multiple  missile  flights  in  sequence 
by  repeating  the  same  deck  organization  shown  above  for  the  multistage 
case.  The  only  difference  is  that  the  variable  IRUN  is  not  specified 
in  BIASIN  after  the  first  run. 

It  should  be  noted  that  values  of  the  input  variables  do  not  carry 
over  from  one  run  to  the  next,  but  revert  to  their  default  values  un- 
less entered  again. 

Block  Data 


In  addition  to  the  input  values  described  above,  there  are  certain 


-53- 


parameters  which  are  set  within  the  program  in  the  block  data  sub- 
program. These  quantities  are  as  follows: 

Degrees  per  radian 

DEG  - 57.2957795131. 

Earth  gravitational  constant 

HMU  - 1.407645  x 1016  ft3/sec1 2 3 4. 

Earth  angular  rate 

OMEG  - 7.29211585  * lO-5  rad/ sec. 

Earth  radius 

AE  - 20,925,696  ft. 

Earth  eccentricity 

EPS  - 0.082. 


OUTPUT  FORMAT 


Title  Page 

On  this  page,  the  program  name  and  a brief  statement  of  the  oper- 
ations performed  by  the  four  major  subroutines  are  printed. 

Input  Parameters  and  Comments 

In  this  section,  the  input  data,  as  well  as  any  initial  conditions 
generated  from  this  input  data,  are  printed  out  in  the  following  order 
at  the  beginning  of  each  stage  of  the  flight. 

1.  The  value  of  the  print  interval,  in  seconds. 

2.  The  values  of  the  flags  ISTG,  IBST,  IMODE,  INPQ,  INPX,  IX, 

IQ,  IS,  together  with  the  interpretation  of  these  numerical 
assignments. 

3.  If  the  covariance  matrix  is  entered  as  diagonal  elements 
only  (INPQ  < 0),  these  diagonal  elements  are  printed  out. 

4.  The  values  of  the  initial  biases  of  the  51  state  variables 
(DXC) , as  selected  from  gaussian  distributions,  or  their  de- 
fault values  of  zero.  These  values  are  only  printed  at  the 
beginning  of  the  first  stage. 
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5.  The  initial  and  terminal  times,  in  seconds,  for  this  stage. 

6.  The  value  of  the  Greenwich  Hour  Angle,  in  degrees,  at  the 
beginning  of  the  stage. 

7.  The  values  of  the  flags  IOBS,  1R,  IH,  IRND,  together  with 
the  interpretation  of  these  numerical  assignments. 

8.  The  initial  position  of  the  missile  and  the  coordinate  sys- 
tem in  which  it  is  expressed.  If  IP  equals  2 or  4,  the  posi- 
tion is  given  in  latitude,  longitude,  and  altitude;  if  IP 
equals  1 or  3,  the  position  is  in  cartesian  coordinates. 

9.  The  initial  missile  velocity  in  the  same  coordinate  system 
as  the  position.  If  IP  equals  2 or  4,  the  velocity  is  given 
in  terms  of  its  path  elevation  angle,  path  azimuth  angle,  and 
velocity  magnitude;  if  IP  equals  1 or  3,  it  is  given  in 
cartesian  coordinates.  These  initial  position  and  velocity 
values  are  printed  out  only  at  the  beginning  of  the  first 
stage. 

10.  The  time  TO  at  which  the  initial  conditions  apply,  measured 
In  seconds,  from  the  start  of  the  run. 

11.  The  maximum  value  of  the  integration  interval  DTP,  in  seconds. 

12.  The  variances  of  state  variables  associated  with  the  IMU. 

13.  The  tabular  acceleration  profile,  giving  the  three  components 
of  missile  acceleration  as  a function  of  time,  if  IBST  does 
not  equal  zero. 

14.  A list  of  which  trackers  are  available  to  make  measurements 
during  this  stage,  where  they  are  based,  and  what  they  can 
measure. 

15.  The  initial  position  and  velocity  of  each  tracker  and  the 
times  at  which  these  conditions  exist. 

16.  The  limits  of  the  intervals  during  which  the  trackers  can 
observe,  the  interval  between  observations,  and  the  maximum 
and  minimum  elevation  angles  of  the  trackers. 

17.  The  measurement  biases  in  range,  range  rate,  elevation  angle, 
and  azimuth  angle  for  each  tracker. 

18.  The  biases  in  latitude,  longitude,  altitude,  and  speed  of 
each  tracker. 
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19.  The  tropospheric  biases  in  range,  range  rate,  and  elevation 
angle  for  an  elevation  angle  of  90  deg. 

20.  The  24  constants  used  to  express  the  standard  deviations  of 
the.  measurements  for  each  tracker. 


Output  Legends 

This  page  is  printed  out  at  the  beginning  of  each  stage  and  gives 
the  format  of  the  data  which  will  be  printed  out  during  this  stage. 

The  format  is  determined  by  the  flags  IQ,  IS,  IX,  IOBS,  IR,  IH,  and 
IPRNT. 


Statistical  Estimation 

This  information  is  printed  out  at  the  end  of  each  print  interval 
and  includes  the  following  items: 

1.  The  time  for  which  these  data  apply,  measured  in  seconds, 
from  the  start  of  the  run. 

2.  The  position  and  velocity  of  the  missile,  expressed  in  the 
coordinate  system,  determined  by  the  value  of  IPRNT. 

3.  The  first  six  elements  of  the  actual  state  vector  DXC  that 
give  the  deviations,  expressed  in  the  coordinate  system 
determined  by  the  flag  IS,  of  the  missile's  position  and 
velocity  from  the  nominal  trajectory. 

4.  The  elements  of  the  estimated  state  vector  DX,  in  the  same 
coordinate  system  as  DXC.  If  no  external  measurements  have 
been  made  since  the  last  print  time,  only  the  first  six 
elements  of  DX  are  printed  for  comparison  with  the  corre- 
sponding DXC  values.  If  measurements  have  been  made,  the 
remaining  45  elements  of  DX  associated  with  the  IMU  param- 
eters are  also  printed  ou*  for  comparison  with  their  actual 
values  listed  in  the  section  on  input  parameters  and  comments. 

5.  The  values  of  the  first  six  elements  of  the  state  variable  X 
that  give  the  current  estimated  position  and  velocity  of  the 
missile  in  the  coordinate  system  specified  by  IS.  If  the 
program  is  being  used  in  the  simulation  mode  (1M0DE=2),  then 
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the  remaining  45  state  variables,  the  estimated  IMU  param- 
eters, are  also  printed  out.  These  values  of  X,  as  well  as 
those  of  DXC  and  DX,  are  only  printed  if  the  flag  IX  is  non- 
zero. 

6.  The  standard  deviations  of  the  estimated  variations  of  51 
state  variables  that  are  the  square  roots  of  the  diagonal 
elements  of  the  state  covariance  matrix. 

7.  The  upper  left  (6  x 6)  of  the  state  covariance  matrix  show- 
ing the  correlation  between  the  components  of  missile  posi- 
tion and  velocity.  The  standard  deviations  and  this  6x6 
matrix  a~e  expressed  in  the  coordinate  system  specified  by 
IS  and  are  only  printed  if  the  *alue  of  the  flag  IQ  is  odd. 

8.  The  remaining  tridiagonal  elements  of  the  state  covariance 
matrix,  showing  the  correlation  between  the  missile  position 
and  velocity  variations,  as  well  as  the  IMU  parameter  varia- 
tions. These  elements  are  printed  if  the  flag  IQ  is  equal 
to  2,  3,  or  7. 

9.  The  eigen  values  and  eigen  vectors  associated  with  the  posi- 
tion covariance  matrix  are  print  in  a six-column  format, 
in  which  the  third  column  is  the  eigen  value  and  the  fourth 
through  sixth  columns  are  the  three  components  of  the  asso- 
ciated eigen  vector.  The  first  column  represents  the  semi- 
axis of  the  50  percent  confidence  ellipsoid,  and  the  second 
column  is  the  semiaxis  of  the  95  percent  confidence  ellipsoid 
for  this  same  vector  direction. 

10.  The  eigen  values  and  eigen  vectors  associated  with  the 
velocity  covariance  matrix  are  printed  out  in  the  same  for- 
mat. The  position  and  velocity  eigen  value  data  are  only 
printed  when  the  flag  IQ  is  greater  than  4. 

11.  The  data  associated  with  observations  by  the  external  trackers 
start  with  a printout  of  the  observation  time,  followed  by 
the  observational  data  of  each  of  the  trackers.  If  a tracker 
makes  no  observation,  the  printout  is  a single  line,  starting 
with  the  tracker  identification  number  and  followed  by  "NONE" 
and  six  zeros.  If  a tracker  does  make  an  observation,  the 
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output  consists  of,  at  most,  four  lines.  The  first  line 
starts  with  the  tracker  identification  number  and  -s  followed 
by  the  latitude,  longitude,  and  altitude  of  the  tracker,  as 
well  as  its  path  elevation  angle,  path  azimuth  angle,  and 
speed.  These  values  are  in  the  same  columns  as  the  six  zeros 
in  the  case  of  the  nonobserving  tracker.  The  second  line 
gives  the  elevation  angle  and  range  of  the  target  relative  to 
the  tracker.  The  third  line  lists  the  nominal  value  of  any 
measurement  made  by  this  tracker  in  columns  between  those 
used  for  the  position  and  velocity  data.  The  order  of  the 
measurements  in  these  columns  is  range,  range  rate,  azimuth 
angle,  elevation  angle,  hyperbolic  range,  and  hyperbolic 
range  rate. 

Alongside  of  each  measurement  value  is  printed  the  total 
number  of  this  type  of  measurement  made  by  this  tracker  up  to 
this  time.  If  IPSMD  or  IPRMD  is  zero,  then  the  minimum  range 
tracker  is  Indicated  in  the  column  for  hyperbolic  range  or 
hyperbolic  range  rate.  If  IPSMD  or  IPRMD  is  unity,  the  -nin- 
imum  range  tracker  is  indicated  by  an  M beside  the  number  of 
hyperbolic  range  or  hyperbolic  range  rate  measurements  for 
that  tracker,  indicating  that  these  are  the  tracker  data  that 
are  propagated  forward  as  a reference. 

The  fourth  line  of  the  printout  for  a given  tracker 
occurs  only  for  the  simulation  mode  of  operation  (IM0DE=2) 
and  gives  the  residual  differences  between  the  measurement 
and  it 8 nominal  value  as  the  reference  trajectory  is  progres- 
sively modified.  These  values  are  printed  in  the  same  columns 
as  the  corresponding  measurements. 

12.  The  elements  of  the  H matrix  are  printed  in  a six-column  for- 
mat with  each  row  corresponding  to  one  measurement,  starting 
with  the  first  measurement  by  the  lowest  numbered  tracker  and 
ending  with  the  last  measurement  by  the  highest  numbered 
tracker.  These  data  are  printed  out  only  if  the  flag  IH  is 


nonzero. 
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13.  The  standard  deviations  of  the  various  measurements  are 
printed  in  a four-column  format  with  one  value  for  each 
measurement  so  that  measurements  1 through  4 are  in  the 
first  row,  5 through  8 in  the  second,  and  so  on — through 
the  final  measurement. 

Expanded  Printout 

If  the  print  interval  PRNTDT  is  larger  than  the  interval  between 
tracker  observations,  the  tracker  observations  and  the  corresponding 
missile  position  and  velocity  data  at  these  intermediate  observation 
times  will  not  be  printed  out.  If  this  information  is  required,  it 
can  be  obtained  by  means  of  the  flags  IPRNT  and  IOBS.  If  IPRNT  is 
negative,  the  missile  position  and  velocity  data  at  each  observation 
time  are  printed  out  in  the  same  form  as  described  in  the  previous 
section  on  statistical  estimation.  By  making  IOBS  negative,  the 
tracking  data  at  each  observation  time  are  printed  out  immediately 
following  the  corresponding  missile  position  data  in  the  same  format 
described  previously. 

At  the  specified  printout  time,  since  the  missile  position  and 
velocity  data  have  already  been  printed,  only  the  values  of  DXC,  DX, 

X,  the  standard  deviations  of  the  estimates,  the  state  covariance 
matrix,  the  eigen  values,  and  the  eigen  vectors  are  printed  out  in 
the  statistical  estimation  format. 

Flight  Termination 

The  final  statistical  estimation  occurs  at  a time  equal  to  either 
TMAX  for  the  final  stage  or  at  the  time  of  impact,  whichever  is  smaller. 
If  the  run  is  terminated  at  TMAX,  the  statistical  estimation  is  in  the 
same  format  as  that  described  above.  However,  if  the  run  is  terminated 
by  Impact  with  earth,  the  following  additional  information  is  printed: 

1.  A statement  indicating  that  the  vehicle  is  down. 

2.  The  semiaxes,  in  feet,  of  the  error  ellipse  on  the  ground  and 
the  orientation  of  its  major  axis  relative  to  an  in-plane 
rotating  coordinate  system. 
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Appendix 

COORDINATE  SYSTEMS 


The  various  coordinate  systems  used  in  the  formulation  of  this 
program  are  described  below,  together  with  the  matrix  transformation 
relating  them. 

EARTH-CENTERED  INERTIAL  (ECI) 

ECI  is  the  basic  coordinate  system,  in  which  all  computations  are 
performed.  Its  origin  is  at  the  center  of  the  earth  with  its  z axis 
colinear  with  the  earth's  spin  axis  and  its  x axis  in  the  direction  of 
the  first  point  of  Aries. 

EARTH-CENTERED  ROTATING  (ECR) 

The  ECR  system  al.  :>  has  its  origin  at  the  earth's  center  and  its 
z axis  along  the  earth's  spin  axis;  its  x axis  is  in  the  plane  of  the 
Greenwich  meridian.  The  relation  between  the  ECR  and  ECI  systems  is 
shown  in  Fig.  1,  where  T is  the  Greenwich  Hour  Angle  between  the  plane 
of  the  Greenwich  meridian  and  the  first  point  of  Aries,  given  by  the 
relation 

r - r0  + «E(t  - tQ)  , (A-i) 


where  Tq  is  the  initial  value  of  T at  some  reference  time  t^,  and  fig 
is  the  spin  angular  rate  of  the  earth. 

The  transformation  of  a vector  from  ECI  to  ECR  coordinates  is 
given  by  the  matrix  equation 


—(ECR) 


r‘eci>  , 


(A-2) 


where  the  transformation  matrix  is  given  by 


ECR 


cos  r sin  r 0 
-sin  T cos  r 0 
0 0 1 
r' 


(A-3) 


T *•“ 


/^B^UffiUNOT  FIIfflSD 


Fig.  1 — Earti  centered  rotating  coordinate  system 


LAUNCH-CENTERED  ROTATING  (LCR) 

The  LCR  system  has  its  origin  at  the  launch  site  of  the  missile 
with  its  z axis  vertical  and  its  x axis  in  the  direction  of  the  launch 
azimuth.  The  relationship  of  this  system  to  the  ECI  system  is  shown 
in  Fig.  2,  where  \Q,  BQ,  and  4*0  are  the  latitude,  longitude,  and  path 
azimuth,  respectively,  at  the  launch  site.  The  transformation  between 
ECI  and  LCR  is  given  by  the  relation 


^(LCR) 


r(eci)  . 
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|*lcr| 


/ -sin  I(i0  sin  (8Q  + T) 
\-cos  sin  Xq  cos  (8 

/ cos  \pQ  sin  (8q  + T) 
y-sin  ijig  sin  XQ  cos  (8 

( cos  XQ  cos  (8q  + r>) 


0 


0 


+ 


+ 


X 

X 


sin  i)(q  cos  (8q  + T) 
-cos  \Jig  sin  Xq  sin  (8 

-cos  cos  (6q  + T) 
-sin  4*0  sin  Xq  sin  (8 

cos  XQ  sin  (8q  + T)) 


0 


0 


.„)< 

.n )( 


C08  C08 


sin  cos 


(sin  XQ) 


xo) 


(A-5) 


Since  F is  a function  of  time  through  Eq.  (A-l),  the  transformation 
matrix  also  changes  with  time  as  the  launch  site  rotates  about  the 
earth's  axis. 

LAUNCH-CENTERED  INERTIAL  (LCI) 

The  LCI  system  also  has  its  origin  at  the  launch  site  with  its 
z axis  in  the  direction  of  the  vertical  at  the  time  of  launch  and 
its  x axis  in  the  direction  of  the  inertial  launch  azimuth.  The 


zeci 


Fig.  2 — Launch  centered  rotating  coordinate  system 


transformation  from  ECI  to  LC1  is  given  by 


r(lci) 


—(ECI) 


(A-6) 


where  the  transformation  matrix  j^jj  has  the  same  form  ao  in 

Eq.  (A-5),  except  that  T is  replaced  by  its  initial  value  Tq.  Thus, 
|alci|  is  independent  of  time  and  the  LCI  coordinate  system  maintains 
a constant  orientation  relative  to  the  ECI  system. 


IN-PLANE  ROTATING  (IPR) 

The  IPR  system  has  its  origin  on  the  missile  or  user  vehicle,  with 
its  z axis  along  the  local  vertical  and  its  x,z  plane  containing  the 
missile  velocity  vector  relative  to  the  rotating  earth.  The  relation 
of  this  system  to  the  ECR  and  ECI  systems  is  shown  in  Fig.  3.  The 


Fig.  3 — In-plane  rotating  coordinate  system 
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transformation  from  ECI  to  IPR  is  given  by 


r(ipr) 


—(ECI) 


(A- 7 ) 


where  the  transformation  matrix  |^jpR|  has  the  same  form  as  that  for 
in  Eq.  (A-5) , with  the  exception  that  the  quantities  Xq,  8q, 
and  iJ/q  are  replaced  by  their  values  at  the  missile's  present  position 
on  its  flight  path. 


IN-PLANE  INERTIAL  (IP I) 

The  IPI  system  differs  from  the  IPR  system  only  in  that  the  x,z 
plane  is  oriented  to  contain  the  velocity  vector  relative  to  inertial 
space.  The  transformation  from  ECI  to  IPI  is  given  by 


R<IPI)  = 


r(eci) 


(A-8) 


where  A has  the  same  form  as  AT1J  except  that  the  value  of  the 

I l I I i*  K I ^ 

path  azimuth  refers  to  the  inertial  velocity  vector. 


PLATFORM  (P) 

This  system  is  associated  with  the  stabilized  platform  on  which 
the  inertial  instruments  are  mounted.  Ideally,  the  platform  coordi- 
nate system  has  the  same  orientation  as  the  LCI  or  LCR  system,  depend- 
ing on  the  particular  guidance  system  used.  Thus,  the  transformation 
from  ECI  to  P is  given  by 

*<p)  - K 


R(ECI) 


(A-9) 


It  should  be  noted  that  the  IPI  system  is  not  actually  an  iner- 
tial reference  system,  since  its  orientation  does  change  relative  to 
inertial  space  although  the  inertial  velocity  remains  in  the  x,z 
plane . 


MISALIGNED  PLATFORM  (P  ) 


■m 


The  actual  platform  has  small  angular  misalignments  relative  to 
dea] 
given  by 


the  idealized  platform  such  that  the  transformation  from  P to  P is 

m 


-(PJ 

R 


m 


R(P> 


(A-l 


where  the  transformation  matrix  is  given  by 


i 

l-K. 


-K, 


m 


K, 


1 

-K, 


(A-l 


and  is  a small  rotation  about  the  yp  axis,  about  the  zp  axis, 
and  K.j  about  the  xp  axis. 


ACCELEROMETER 

This  system  is  associated  with  the  ith  accelerometer  with  its 
x axis  along  the  input  axis,  y along  the  cross  axis,  and  z along  the 
normal  axis.  Its  orientation  relative  to  the  P system  is  shown  in 
Fig.  4,  and  involves  a rotation  through  an  angle  ot^  about  the  y axis, 
followed  by  a rotation  B^  about  the  z axis,  and  finally  a rotation  of 
about  the  x axis.  The  resulting  transformation  from  P to  a^  is 
given  by 


R 


-67- 


zp 


Fig,  4 — Accelerometer  coordinate  system 


where  the  transformation  matrix  is  given  by 


COB  cos 


h ) ( sln  ei  ) 


(-sin  cos 


o 


/ Bln  sin  Y*  \/  co8  B c08  Y \ / cos  al  8in  Yi  \ 

\-cos  at  sin  flt  cos  Y^V  7 W"  ai  8in  61  cos  Yi  / 

/ 8ln  °l  C08  Yi  \ /-cos  6t  sin  y\  ( co8  ai  C08  yi  ) 

y+cos  sin  sin  Y±/  \ /\-sln  sin  61  sin  Y±  / 


(A- 14) 


GYROSCOPE  (gi) 

The  gt  system  is  associated  with  the  ith  gyroscope  and  has  its 
x axis  along  the  input  axis,  y along  the  output  axis,  and  z along  the 
spin  axis.  As  in  the  case  of  the  accelerometer  axes,  the  orientation 
of  relative  to  the  platform  axes  involves  rotations  of  a^,  8^,  and 


i 
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tracker  longitude.  The  transformation  from  ECR  to  T is  given  hy 


(T.) 


(ECR) 


(A-16) 


where  the  transformation  matrix  has  the  form 
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(A-17) 


OTHER  TRANSFORMATIONS 

The  transformation  matrix  between  any  two  of  the  coordinate  sys- 
tems defined  above  can  be  obtained  as  the  products  of  the  matrices 
defined  above  and  their  inverses. 


STATE  VECTOR  TRANSFORMATION  (A) 

The  matrix  [A]  indicated  in  Eq.  (39)  represents  the  transforma- 
tion of  the  position  and  velocity  elements  of  the  state  vector  from 
ECI  coordinates  to  tracker  coordinates.  The  derivation  of  this  matrix 
is  as  follows.  If  the  position  vector  in  tracker  coordinates  is  given 
as 


(A-18) 


and  the  corresponding  position  vector  in  ECI  coordinates  is 


if(ECI) 


(A-19) 
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then  the  two  are  related  by  the  equation 


_cr.) 

R 


— (EC1) 


1 o 

1 

J 

* where  the  matrix  1 A 

1 ° 

is 

given  by 

A 1 = |A_ 

X 

A 

-1  N 

'! 

ECR 

(A-20) 


(A-21) 


The  corresponding  relations  between  the  velocity  vectors  in  the  two 
systems  are  obtained  as  follows. 

If  the  velocity  vector  in  tracker  coordinates  is  given  by 


(A-22) 


and  the  corresponding  vector  in  EC1  coordinates  is 


r(ECI)  = 


(A-23) 


then  the  two  are  related  by  the  equation 


(A-24) 


(rrj) 

where  wv  represents  the  angular  velocity  of  the  tracker  coordinate 
system  relative  to  inertial  space  in  ECI  coordinates.  This  angular 
velocity  is  the  result  of  the  rotation  of  the  earth  and  the  motion  of 
the  tracker  relative  to  the  rotating  earth.  In  ECR  coordinates,  this 
angular  velocity  is  given  by 
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-(ECR) 


-ecr 

WE 


(A-25) 


/ irr’D  j 

where  R^,  is  the  vector  position  of  the  tracker  in  ECR  coordinates 
given  by  the  expression 


-(ECR)  . 


(A-26) 


—(VCR) 

and  is  the  velocity  of  the  tracker  relative  to  the  earth  in 

ECR  coordinates  expressed  as 


~(ECR) 

*T 


(A-27 ) 


Equation  (A-25)  can  be  rewritten  in  matrix  form  as 


-(ECR) 

U)  = 
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0 ~ZT  yT 

*T 

0 

ZT  ° -XT 

X 

yT 

*1 

WE 

*T  0 

£t 

(yTzT  - zTyT)/RT 


(zt*t  - Vi)/rt 

" yT*T^RT  + UE 


(A-28) 


This  angular  rate  can  then  be  transformed  into  EC1  coordinates  as 


r-^  ( 
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